diff --git a/src/Solution/ParticleTracker/Cell.f90 b/src/Solution/ParticleTracker/Cell.f90 index 36d95c207eb..a1e2b8e7a1b 100644 --- a/src/Solution/ParticleTracker/Cell.f90 +++ b/src/Solution/ParticleTracker/Cell.f90 @@ -8,7 +8,7 @@ module CellModule !> @brief A grid cell. Contains a cell-definition (composition over inheritance) type, abstract :: CellType character(len=40), pointer :: type ! tracking domain type - type(CellDefnType), pointer :: cellDefn => null() ! pointer to cell definition, todo rename to defn + type(CellDefnType), pointer :: defn => null() ! cell definition contains procedure(destroy), deferred :: destroy ! destructor for the cell end type CellType diff --git a/src/Solution/ParticleTracker/CellPoly.f90 b/src/Solution/ParticleTracker/CellPoly.f90 index 407057efc5e..a0aade23a2b 100644 --- a/src/Solution/ParticleTracker/CellPoly.f90 +++ b/src/Solution/ParticleTracker/CellPoly.f90 @@ -19,7 +19,7 @@ module CellPolyModule subroutine create_cellPoly(cellPoly) type(CellPolyType), pointer :: cellPoly allocate (cellPoly) - allocate (cellPoly%cellDefn) + allocate (cellPoly%defn) allocate (cellPoly%type) cellPoly%type = 'CellPoly' end subroutine create_cellPoly @@ -27,7 +27,7 @@ end subroutine create_cellPoly !> @brief Destroy the polygonal cell subroutine destroy_cellPoly(this) class(CellPolyType), intent(inout) :: this - deallocate (this%cellDefn) + deallocate (this%defn) deallocate (this%type) end subroutine destroy_cellPoly diff --git a/src/Solution/ParticleTracker/CellRect.f90 b/src/Solution/ParticleTracker/CellRect.f90 index bd39f6ed363..e234f84a28a 100644 --- a/src/Solution/ParticleTracker/CellRect.f90 +++ b/src/Solution/ParticleTracker/CellRect.f90 @@ -35,7 +35,7 @@ module CellRectModule subroutine create_cellRect(cellRect) type(CellRectType), pointer :: cellRect allocate (cellRect) - allocate (cellRect%cellDefn) + allocate (cellRect%defn) allocate (cellRect%type) cellRect%type = 'CellRect' end subroutine create_cellRect @@ -43,7 +43,7 @@ end subroutine create_cellRect !> @brief Destructor for a rectangular cell subroutine destroy_cellRect(this) class(CellRectType), intent(inout) :: this - deallocate (this%cellDefn) + deallocate (this%defn) deallocate (this%type) end subroutine destroy_cellRect diff --git a/src/Solution/ParticleTracker/CellRectQuad.f90 b/src/Solution/ParticleTracker/CellRectQuad.f90 index a58e48f7f50..5d4e0d70db4 100644 --- a/src/Solution/ParticleTracker/CellRectQuad.f90 +++ b/src/Solution/ParticleTracker/CellRectQuad.f90 @@ -41,7 +41,7 @@ module CellRectQuadModule subroutine create_cellRectQuad(cellRectQuad) type(CellRectQuadType), pointer :: cellRectQuad allocate (cellRectQuad) - allocate (cellRectQuad%cellDefn) + allocate (cellRectQuad%defn) allocate (cellRectQuad%irectvert(5)) allocate (cellRectQuad%ipv4irv(2, 4)) allocate (cellRectQuad%rectflow(2, 4)) @@ -52,7 +52,7 @@ end subroutine create_cellRectQuad !> @brief Destroy a rectangular-quad cell subroutine destroy_cellRectQuad(this) class(CellRectQuadType), intent(inout) :: this - deallocate (this%cellDefn) + deallocate (this%defn) deallocate (this%irectvert) deallocate (this%type) end subroutine destroy_cellRectQuad @@ -64,7 +64,7 @@ subroutine init_from(this, cellDefn) type(CellDefnType), pointer :: cellDefn ! -- Set pointer to cell definition - this%cellDefn => cellDefn + this%defn => cellDefn ! -- Load the "rectangle vertices" for the cell call this%load_irectvert() @@ -81,21 +81,21 @@ subroutine load_irectvert(this) ! -- local integer :: npolyverts, n, m - npolyverts = this%cellDefn%get_npolyverts() + npolyverts = this%defn%get_npolyverts() n = 0 do m = 1, npolyverts - if (.not. this%cellDefn%get_ispv180(m)) then + if (.not. this%defn%get_ispv180(m)) then n = n + 1 this%irectvert(n) = m this%ipv4irv(1, n) = m - this%rectflow(1, n) = this%cellDefn%get_faceflow(m) + this%rectflow(1, n) = this%defn%get_faceflow(m) this%ipv4irv(2, n) = 0 this%rectflow(2, n) = 0d0 else if (n .ne. 0) then this%ipv4irv(2, n) = m - this%rectflow(2, n) = this%cellDefn%get_faceflow(m) + this%rectflow(2, n) = this%defn%get_faceflow(m) end if end if end do @@ -130,19 +130,19 @@ function get_irectvertSW(this) result(irv1) irv1 = irvnxt(irv4) ipv4 = this%irectvert(irv4) ipv1 = this%irectvert(irv1) - x4 = this%cellDefn%polyvert(1, ipv4) - y4 = this%cellDefn%polyvert(2, ipv4) - x1 = this%cellDefn%polyvert(1, ipv1) - y1 = this%cellDefn%polyvert(2, ipv1) + x4 = this%defn%polyvert(1, ipv4) + y4 = this%defn%polyvert(2, ipv4) + x1 = this%defn%polyvert(1, ipv1) + y1 = this%defn%polyvert(2, ipv1) if (x1 .lt. x4) then irv2 = irvnxt(irv1) ipv2 = this%irectvert(irv2) - x2 = this%cellDefn%polyvert(1, ipv2) + x2 = this%defn%polyvert(1, ipv2) if (x2 .le. x1) return else if (y1 .ge. y4) then irv2 = irvnxt(irv1) ipv2 = this%irectvert(irv2) - y2 = this%cellDefn%polyvert(2, ipv2) + y2 = this%defn%polyvert(2, ipv2) if (y2 .gt. y1) return end if end do @@ -171,26 +171,26 @@ subroutine get_rectDimensionsRotation(this, irv1, xOrigin, yOrigin, zOrigin, & ! -- Get model coordinates at irv1, irv2, and irv4 ipv1 = this%irectvert(irv1) - x1 = this%cellDefn%polyvert(1, ipv1) - y1 = this%cellDefn%polyvert(2, ipv1) + x1 = this%defn%polyvert(1, ipv1) + y1 = this%defn%polyvert(2, ipv1) ipv2 = this%irectvert(irv2) - x2 = this%cellDefn%polyvert(1, ipv2) - y2 = this%cellDefn%polyvert(2, ipv2) + x2 = this%defn%polyvert(1, ipv2) + y2 = this%defn%polyvert(2, ipv2) ipv4 = this%irectvert(irv4) - x4 = this%cellDefn%polyvert(1, ipv4) - y4 = this%cellDefn%polyvert(2, ipv4) + x4 = this%defn%polyvert(1, ipv4) + y4 = this%defn%polyvert(2, ipv4) ! -- Compute rectangle dimensions xOrigin = x1 yOrigin = y1 - zOrigin = this%cellDefn%bot + zOrigin = this%defn%bot dx2 = x2 - xOrigin dy2 = y2 - yOrigin dx4 = x4 - xOrigin dy4 = y4 - yOrigin dx = dsqrt(dx4 * dx4 + dy4 * dy4) dy = dsqrt(dx2 * dx2 + dy2 * dy2) - dz = this%cellDefn%top - zOrigin ! kluge note: need to account for partial saturation + dz = this%defn%top - zOrigin ! kluge note: need to account for partial saturation ! -- Compute sine and cosine of rotation angle (angle between "southern" ! -- rectangle side irv1-irv4 and the model x axis) diff --git a/src/Solution/ParticleTracker/CellUtil.f90 b/src/Solution/ParticleTracker/CellUtil.f90 index adfa4676919..5e0d62f2770 100644 --- a/src/Solution/ParticleTracker/CellUtil.f90 +++ b/src/Solution/ParticleTracker/CellUtil.f90 @@ -29,7 +29,7 @@ subroutine CellPolyToCellRect(cellPoly, cellRect, istatus) double precision :: factor, term call create_cellRect(cellRect) - cellDefn => cellPoly%cellDefn + cellDefn => cellPoly%defn ! -- kluge note: no check whether conversion is possible; assumes it is ! -- Translate and rotate the rectangular cell into local coordinates @@ -85,7 +85,7 @@ subroutine CellPolyToCellRect(cellPoly, cellRect, istatus) dz = cellDefn%top - zOrigin ! todo: need to account for partial saturation sinrot = dy4 / dx cosrot = dx4 / dx - cellRect%cellDefn = cellPoly%cellDefn + cellRect%defn = cellPoly%defn cellRect%dx = dx cellRect%dy = dy cellRect%dz = dz @@ -129,7 +129,7 @@ subroutine CellPolyToCellRectQuad(cellPoly, cellRectQuad, istatus) double precision :: qhalf, qdisttopbot, q1, q2, q4 call create_cellRectQuad(cellRectQuad) - call cellRectQuad%init_from(cellPoly%cellDefn) + call cellRectQuad%init_from(cellPoly%defn) ! kluge note: no check whether conversion is possible; assumes it is ! -- Translate and rotate the rect-quad cell into local coordinates with ! -- x varying from 0 to dx and y varying from 0 to dy. Choose the "south- @@ -159,9 +159,9 @@ subroutine CellPolyToCellRectQuad(cellPoly, cellRectQuad, istatus) cellRectQuad%qextl1(isc) = cellRectQuad%get_rectflow(2, irv) end if end do - qdisttopbot = 2.5d-1 * (cellRectQuad%cellDefn%get_distflow() & - + cellRectQuad%cellDefn%get_botflow() & - + cellRectQuad%cellDefn%get_topflow()) + qdisttopbot = 2.5d-1 * (cellRectQuad%defn%get_distflow() & + + cellRectQuad%defn%get_botflow() & + + cellRectQuad%defn%get_topflow()) q1 = qdisttopbot + cellRectQuad%qextl1(1) + cellRectQuad%qextl2(1) q2 = qdisttopbot + cellRectQuad%qextl1(2) + cellRectQuad%qextl2(2) q4 = qdisttopbot + cellRectQuad%qextl1(4) + cellRectQuad%qextl2(4) diff --git a/src/Solution/ParticleTracker/MethodCellPollock.f90 b/src/Solution/ParticleTracker/MethodCellPollock.f90 index 84dc90ee7d3..0a3ad370aa2 100644 --- a/src/Solution/ParticleTracker/MethodCellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollock.f90 @@ -137,7 +137,7 @@ subroutine apply_mCP(this, particle, tmax) ! -- Update particle state, checking whether any reporting or ! -- termination conditions apply - call this%update(particle, this%cellRect%cellDefn) + call this%update(particle, this%cellRect%defn) ! -- Return early if particle is done advancing if (.not. particle%advancing) return @@ -146,8 +146,8 @@ subroutine apply_mCP(this, particle, tmax) ! -- represent a water table above the cell bottom), pass the particle ! -- vertically and instantaneously to the cell top elevation and save ! -- the particle state to output file(s). - if (particle%z > this%cellRect%cellDefn%top) then - particle%z = this%cellRect%cellDefn%top + if (particle%z > this%cellRect%defn%top) then + particle%z = this%cellRect%defn%top call this%trackctl%save(particle, kper=kper, & kstp=kstp, reason=1) ! reason=1: cell transition end if diff --git a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 index 8342058c3da..aa9c5c50d67 100644 --- a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 @@ -94,7 +94,7 @@ subroutine pass_mCPQ(this, particle) exitFace = particle%iTrackingDomainBoundary(3) isc = particle%iTrackingDomain(3) - npolyverts = this%cellRectQuad%cellDefn%npolyverts + npolyverts = this%cellRectQuad%defn%npolyverts select case (exitFace) ! kluge note: exitFace uses Dave's iface convention case (0) @@ -221,7 +221,7 @@ subroutine apply_mCPQ(this, particle, tmax) ! -- Update particle state, checking whether any reporting or ! -- termination conditions apply - call this%update(particle, this%cellRectQuad%cellDefn) + call this%update(particle, this%cellRectQuad%defn) ! -- Return early if particle is done advancing if (.not. particle%advancing) return @@ -230,8 +230,8 @@ subroutine apply_mCPQ(this, particle, tmax) ! -- represent a water table above the cell bottom), pass the particle ! -- vertically and instantaneously to the cell top elevation and save ! -- the particle state to output file(s). - if (particle%z > this%cellRectQuad%cellDefn%top) then - particle%z = this%cellRectQuad%cellDefn%top + if (particle%z > this%cellRectQuad%defn%top) then + particle%z = this%cellRectQuad%defn%top call this%trackctl%save(particle, kper=kper, & kstp=kstp, reason=1) ! reason=1: cell transition end if @@ -274,9 +274,9 @@ subroutine load_subcell(this, particle, levelNext, subcellRect) double precision :: qextl1, qextl2, qintl1, qintl2 double precision :: factor, term - factor = DONE / this%cellRectQuad%cellDefn%retfactor - factor = factor / this%cellRectQuad%cellDefn%porosity - npolyverts = this%cellRectQuad%cellDefn%npolyverts + factor = DONE / this%cellRectQuad%defn%retfactor + factor = factor / this%cellRectQuad%defn%porosity + npolyverts = this%cellRectQuad%defn%npolyverts isc = particle%iTrackingDomain(3) ! -- Subcells 1, 2, 3, and 4 are Pollock's subcells A, B, C, and D, @@ -324,8 +324,8 @@ subroutine load_subcell(this, particle, levelNext, subcellRect) end if dx = 5d-1 * dx dy = 5d-1 * dy - dz = this%cellRectQuad%cellDefn%top - & - this%cellRectQuad%cellDefn%bot ! kluge note: need to account for partial saturation + dz = this%cellRectQuad%defn%top - & + this%cellRectQuad%defn%bot ! kluge note: need to account for partial saturation areax = dy * dz areay = dx * dz areaz = dx * dy @@ -382,8 +382,8 @@ subroutine load_subcell(this, particle, levelNext, subcellRect) m1 = npolyverts + 2 m2 = m1 + 1 term = factor / areaz - subcellRect%vz1 = 2.5d-1 * this%cellRectQuad%cellDefn%faceflow(m1) * term - subcellRect%vz2 = -2.5d-1 * this%cellRectQuad%cellDefn%faceflow(m2) * term + subcellRect%vz1 = 2.5d-1 * this%cellRectQuad%defn%faceflow(m1) * term + subcellRect%vz2 = -2.5d-1 * this%cellRectQuad%defn%faceflow(m2) * term end subroutine load_subcell diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index 1fd7ee6ff2b..a32f86da94a 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -93,7 +93,7 @@ subroutine pass_mCT(this, particle) exitFace = particle%iTrackingDomainBoundary(3) isc = particle%iTrackingDomain(3) - npolyverts = this%cellPoly%cellDefn%npolyverts + npolyverts = this%cellPoly%defn%npolyverts select case (exitFace) case (0) @@ -150,7 +150,7 @@ subroutine apply_mCT(this, particle, tmax) ! -- Update particle state, checking whether any reporting or ! -- termination conditions apply - call this%update(particle, this%cellPoly%cellDefn) + call this%update(particle, this%cellPoly%defn) ! -- Return early if particle is done advancing if (.not. particle%advancing) return @@ -159,36 +159,36 @@ subroutine apply_mCT(this, particle, tmax) ! -- represent a water table above the cell bottom), pass the particle ! -- vertically and instantaneously to the cell top elevation and save ! -- the particle state to output file(s). - if (particle%z > this%cellPoly%cellDefn%top) then - particle%z = this%cellPoly%cellDefn%top + if (particle%z > this%cellPoly%defn%top) then + particle%z = this%cellPoly%defn%top call this%trackctl%save(particle, kper=kper, & kstp=kstp, reason=1) ! reason=1: cell transition end if - npolyverts = this%cellPoly%cellDefn%npolyverts + npolyverts = this%cellPoly%defn%npolyverts xsum = DZERO ysum = DZERO vxsum = DZERO vysum = DZERO area = DZERO - this%ztop = this%cellPoly%cellDefn%top - this%zbot = this%cellPoly%cellDefn%bot + this%ztop = this%cellPoly%defn%top + this%zbot = this%cellPoly%defn%bot this%dz = this%ztop - this%zbot do iv = 1, npolyverts ivp1 = iv + 1 if (ivp1 .gt. npolyverts) ivp1 = 1 ivm1 = iv - 1 if (ivm1 .lt. 1) ivm1 = npolyverts - x0 = this%cellPoly%cellDefn%polyvert(1, iv) - y0 = this%cellPoly%cellDefn%polyvert(2, iv) - x2 = this%cellPoly%cellDefn%polyvert(1, ivp1) - y2 = this%cellPoly%cellDefn%polyvert(2, ivp1) - x1 = this%cellPoly%cellDefn%polyvert(1, ivm1) - y1 = this%cellPoly%cellDefn%polyvert(2, ivm1) - term = DONE / (this%cellPoly%cellDefn%porosity * this%dz) - flow0 = this%cellPoly%cellDefn%faceflow(iv) * term - flow1 = this%cellPoly%cellDefn%faceflow(ivm1) * term + x0 = this%cellPoly%defn%polyvert(1, iv) + y0 = this%cellPoly%defn%polyvert(2, iv) + x2 = this%cellPoly%defn%polyvert(1, ivp1) + y2 = this%cellPoly%defn%polyvert(2, ivp1) + x1 = this%cellPoly%defn%polyvert(1, ivm1) + y1 = this%cellPoly%defn%polyvert(2, ivm1) + term = DONE / (this%cellPoly%defn%porosity * this%dz) + flow0 = this%cellPoly%defn%faceflow(iv) * term + flow1 = this%cellPoly%defn%faceflow(ivm1) * term d01x = x1 - x0 ! kluge note: do this more efficiently, not recomputing things so much??? d01y = y1 - y0 d02x = x2 - x0 @@ -200,7 +200,7 @@ subroutine apply_mCT(this, particle, tmax) ! v0x = -velmult*oodet*(d02x*flow1 + d01x*flow0) ! v0y = -velmult*oodet*(d02y*flow1 + d01y*flow0) ! det = d01y * d02x - d02y * d01x - retfactor = this%cellPoly%cellDefn%retfactor + retfactor = this%cellPoly%defn%retfactor ! kluge note: can det ever be zero, like maybe for a 180-deg vertex??? ! term = velfactor/det ! kluge note: can det ever be zero, like maybe for a 180-deg vertex??? @@ -220,9 +220,9 @@ subroutine apply_mCT(this, particle, tmax) area = area + x0 * y1 - x1 * y0 end do area = area * DHALF - term = DONE / (retfactor * this%cellPoly%cellDefn%porosity * area) - this%vzbot = this%cellPoly%cellDefn%faceflow(npolyverts + 2) * term - this%vztop = -this%cellPoly%cellDefn%faceflow(npolyverts + 3) * term + term = DONE / (retfactor * this%cellPoly%defn%porosity * area) + this%vzbot = this%cellPoly%defn%faceflow(npolyverts + 2) * term + this%vztop = -this%cellPoly%defn%faceflow(npolyverts + 3) * term this%xctr = xsum / dble(npolyverts) this%yctr = ysum / dble(npolyverts) this%vxctr = vxsum / dble(npolyverts) @@ -249,10 +249,10 @@ subroutine load_subcell(this, particle, levelNext, subcellTri) double precision :: di2, d02, d12, di1, d01, alphai, betai double precision :: betatol - ic = this%cellPoly%cellDefn%icell + ic = this%cellPoly%defn%icell subcellTri%icell = ic isc = particle%iTrackingDomain(3) - npolyverts = this%cellPoly%cellDefn%npolyverts + npolyverts = this%cellPoly%defn%npolyverts ! -- Find subcell if not known ! kluge note: from "find_init_triangle" if (isc .le. 0) then diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90 index f174cb4bf63..ac1603ca894 100644 --- a/src/Solution/ParticleTracker/MethodDis.f90 +++ b/src/Solution/ParticleTracker/MethodDis.f90 @@ -102,12 +102,12 @@ subroutine loadsub_mGD(this, particle, levelNext, submethod) double precision :: factor, term ic = particle%iTrackingDomain(levelNext) ! kluge note: is cell number always known coming in? - call this%load_cellDefn(ic, this%cellRect%cellDefn) + call this%load_cellDefn(ic, this%cellRect%defn) ! -- If cell is active but dry, select and initialize ! -- pass-to-bottom method and set cell method pointer if (this%fmi%ibdgwfsat0(ic) == 0) then ! kluge note: use cellDefn%sat == DZERO here instead? - call methodCellPassToBot%init(particle, this%cellRect%cellDefn, & + call methodCellPassToBot%init(particle, this%cellRect%defn, & this%trackctl) submethod => methodCellPassToBot else @@ -119,30 +119,30 @@ subroutine loadsub_mGD(this, particle, levelNext, submethod) dx = dis%delr(jcol) dy = dis%delc(irow) end select - dz = this%cellRect%cellDefn%top - this%cellRect%cellDefn%bot + dz = this%cellRect%defn%top - this%cellRect%defn%bot this%cellRect%dx = dx this%cellRect%dy = dy this%cellRect%dz = dz this%cellRect%sinrot = DZERO this%cellRect%cosrot = DONE - this%cellRect%xOrigin = this%cellRect%cellDefn%polyvert(1, 1) ! kluge note: could avoid using polyvert here - this%cellRect%yOrigin = this%cellRect%cellDefn%polyvert(2, 1) - this%cellRect%zOrigin = this%cellRect%cellDefn%bot + this%cellRect%xOrigin = this%cellRect%defn%polyvert(1, 1) ! kluge note: could avoid using polyvert here + this%cellRect%yOrigin = this%cellRect%defn%polyvert(2, 1) + this%cellRect%zOrigin = this%cellRect%defn%bot this%cellRect%ipvOrigin = 1 areax = dy * dz areay = dx * dz areaz = dx * dy - factor = DONE / this%cellRect%cellDefn%retfactor - factor = factor / this%cellRect%cellDefn%porosity + factor = DONE / this%cellRect%defn%retfactor + factor = factor / this%cellRect%defn%porosity term = factor / areax - this%cellRect%vx1 = this%cellRect%cellDefn%faceflow(1) * term - this%cellRect%vx2 = -this%cellRect%cellDefn%faceflow(3) * term + this%cellRect%vx1 = this%cellRect%defn%faceflow(1) * term + this%cellRect%vx2 = -this%cellRect%defn%faceflow(3) * term term = factor / areay - this%cellRect%vy1 = this%cellRect%cellDefn%faceflow(4) * term - this%cellRect%vy2 = -this%cellRect%cellDefn%faceflow(2) * term + this%cellRect%vy1 = this%cellRect%defn%faceflow(4) * term + this%cellRect%vy2 = -this%cellRect%defn%faceflow(2) * term term = factor / areaz - this%cellRect%vz1 = this%cellRect%cellDefn%faceflow(6) * term - this%cellRect%vz2 = -this%cellRect%cellDefn%faceflow(7) * term + this%cellRect%vz1 = this%cellRect%defn%faceflow(6) * term + this%cellRect%vz2 = -this%cellRect%defn%faceflow(7) * term ! -- Select and initialize Pollock's method and set method pointer call methodCellPollock%init(particle, this%cellRect, this%trackctl) @@ -166,7 +166,7 @@ subroutine pass_mGD(this, particle) inface = particle%iTrackingDomainBoundary(2) z = particle%z - inbr = this%cellRect%cellDefn%facenbr(inface) + inbr = this%cellRect%defn%facenbr(inface) if (inbr .eq. 0) then ! -- Exterior face; no neighbor to map to ! particle%iTrackingDomain(1) = 0 @@ -179,7 +179,7 @@ subroutine pass_mGD(this, particle) kstp=kstp, reason=3) ! reason=3: termination ! particle%iTrackingDomainBoundary(2) = -1 else - idiag = this%fmi%dis%con%ia(this%cellRect%cellDefn%icell) + idiag = this%fmi%dis%con%ia(this%cellRect%defn%icell) ipos = idiag + inbr ic = this%fmi%dis%con%ja(ipos) ! kluge note: use PRT model's DIS instead of fmi's??? particle%iTrackingDomain(2) = ic @@ -210,8 +210,8 @@ subroutine pass_mGD(this, particle) particle%iTrackingDomainBoundary(2) = inface if (inface < 5) then ! -- Map z between cells - topfrom = this%cellRect%cellDefn%top - botfrom = this%cellRect%cellDefn%bot + topfrom = this%cellRect%defn%top + botfrom = this%cellRect%defn%bot zrel = (z - botfrom) / (topfrom - botfrom) top = this%fmi%dis%top(ic) ! kluge note: use PRT model's DIS instead of fmi's??? bot = this%fmi%dis%bot(ic) @@ -259,101 +259,101 @@ function get_top(this, iatop) result(top) end function get_top !> @brief Loads cell definition from the grid - subroutine load_cellDefn(this, ic, cellDefn) + subroutine load_cellDefn(this, ic, defn) ! -- dummy class(MethodDisType), intent(inout) :: this integer, intent(in) :: ic - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- Load basic cell definition components - call this%load_cellDefn_basic(ic, cellDefn) + call this%load_cellDefn_basic(ic, defn) ! -- Load polygon vertices - call this%load_cellDefn_polyverts(cellDefn) + call this%load_cellDefn_polyverts(defn) ! -- Load face neighbors - call this%load_cellDefn_facenbr(cellDefn) + call this%load_cellDefn_facenbr(defn) ! -- Load 180-degree indicator - call this%load_cellDefn_ispv180(cellDefn) + call this%load_cellDefn_ispv180(defn) ! -- Load flows (assumes face neighbors already loaded) - call this%load_cellDefn_flows(cellDefn) + call this%load_cellDefn_flows(defn) end subroutine load_cellDefn !> @brief Loads basic cell definition components from the grid. !! These include cell index, vertex count, (ia)top and botm, !! porosity, retfactor, and izone. - subroutine load_cellDefn_basic(this, ic, cellDefn) + subroutine load_cellDefn_basic(this, ic, defn) ! -- dummy class(MethodDisType), intent(inout) :: this integer, intent(in) :: ic - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- local integer :: iatop real(DP) :: top, bot, sat ! -- cell index - cellDefn%icell = ic + defn%icell = ic ! -- number of polygon vertices - cellDefn%npolyverts = 4 ! rectangular cell always has 4 vertices + defn%npolyverts = 4 ! rectangular cell always has 4 vertices ! -- iatop, top, and bot iatop = get_iatop(this%fmi%dis%get_ncpl(), & this%fmi%dis%get_nodeuser(ic)) - cellDefn%iatop = iatop + defn%iatop = iatop top = this%fmi%dis%top(ic) bot = this%fmi%dis%bot(ic) sat = this%fmi%gwfsat(ic) top = bot + sat * (top - bot) - cellDefn%top = top - cellDefn%bot = bot - cellDefn%sat = sat + defn%top = top + defn%bot = bot + defn%sat = sat ! -- porosity, retfactor, and izone - cellDefn%porosity = this%porosity(ic) - cellDefn%retfactor = this%retfactor(ic) - cellDefn%izone = this%izone(ic) + defn%porosity = this%porosity(ic) + defn%retfactor = this%retfactor(ic) + defn%izone = this%izone(ic) end subroutine load_cellDefn_basic !> @brief Load polygon vertices into cell definition from the grid !! Assumes cell index is already loaded into cell definition. - subroutine load_cellDefn_polyverts(this, cellDefn) + subroutine load_cellDefn_polyverts(this, defn) ! -- dummy class(MethodDisType), intent(inout) :: this - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- Load polygon vertices call this%fmi%dis%get_polyverts( & - cellDefn%icell, & - cellDefn%polyvert, & + defn%icell, & + defn%polyvert, & closed=.true.) end subroutine load_cellDefn_polyverts !> @brief Loads face neighbors to cell definition from the grid. !! Assumes cell index and number of vertices are already loaded. - subroutine load_cellDefn_facenbr(this, cellDefn) + subroutine load_cellDefn_facenbr(this, defn) use GwfDisModule ! -- dummy class(MethodDisType), intent(inout) :: this - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- local integer :: ic1, ic2, icu1, icu2, j1, iloc, ipos integer :: irow1, irow2, jcol1, jcol2, klay1, klay2 integer :: iedgeface ! -- Allocate facenbr array - call ExpandArray(cellDefn%facenbr, cellDefn%npolyverts + 3) + call ExpandArray(defn%facenbr, defn%npolyverts + 3) ! -- Load face neighbors select type (dis => this%fmi%dis) ! kluge type guard type is (GwfDisType) ! kluge - cellDefn%facenbr = 0 - ic1 = cellDefn%icell + defn%facenbr = 0 + ic1 = defn%icell icu1 = dis%get_nodeuser(ic1) call get_ijk(icu1, dis%nrow, dis%ncol, dis%nlay, irow1, jcol1, klay1) call get_jk(icu1, dis%get_ncpl(), dis%nlay, j1, klay1) @@ -378,18 +378,18 @@ subroutine load_cellDefn_facenbr(this, cellDefn) ! Neighbor to the W iedgeface = 1 end if - cellDefn%facenbr(iedgeface) = int(iloc, 1) + defn%facenbr(iedgeface) = int(iloc, 1) else if (klay2 > klay1) then ! -- Bottom face neighbor - cellDefn%facenbr(cellDefn%npolyverts + 2) = int(iloc, 1) + defn%facenbr(defn%npolyverts + 2) = int(iloc, 1) else ! -- Top face neighbor - cellDefn%facenbr(cellDefn%npolyverts + 3) = int(iloc, 1) + defn%facenbr(defn%npolyverts + 3) = int(iloc, 1) end if end do ! -- List of edge (polygon) faces wraps around ! todo: why need to wrap around? no analog to "closing" a polygon? - cellDefn%facenbr(cellDefn%npolyverts + 1) = cellDefn%facenbr(1) + defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1) end select end subroutine load_cellDefn_facenbr @@ -397,78 +397,78 @@ end subroutine load_cellDefn_facenbr !> @brief Load flows into the cell definition. !! These include face flows and net distributed flows. !! Assumes cell index and number of vertices are already loaded. - subroutine load_cellDefn_flows(this, cellDefn) + subroutine load_cellDefn_flows(this, defn) ! -- dummy class(MethodDisType), intent(inout) :: this - type(CellDefnType), intent(inout) :: cellDefn + type(CellDefnType), intent(inout) :: defn ! -- local integer :: ic, m, n, npolyverts - ic = cellDefn%icell - npolyverts = cellDefn%npolyverts + ic = defn%icell + npolyverts = defn%npolyverts ! -- allocate faceflow array - call ExpandArray(cellDefn%faceflow, npolyverts + 3) + call ExpandArray(defn%faceflow, npolyverts + 3) ! -- Load face flows. Note that the faceflow array ! -- does not get reallocated if it is already allocated ! -- to a size greater than or equal to npolyverts+3. - cellDefn%faceflow = 0d0 ! kluge note: eventually use DZERO for 0d0 throughout + defn%faceflow = 0d0 ! kluge note: eventually use DZERO for 0d0 throughout ! -- As with polygon nbrs, polygon face flows wrap around for ! -- convenience at position npolyverts+1, and bot and top flows ! -- are tacked on the end of the list do m = 1, npolyverts + 3 - n = cellDefn%facenbr(m) + n = defn%facenbr(m) if (n > 0) & - cellDefn%faceflow(m) = this%fmi%gwfflowja(this%fmi%dis%con%ia(ic) + n) - ! if (cellDefn%faceflow(m) < 0d0) cellDefn%inoexitface = 0 + defn%faceflow(m) = this%fmi%gwfflowja(this%fmi%dis%con%ia(ic) + n) + ! if (cellDefn%faceflow(m) < 0d0) defn%inoexitface = 0 end do ! -- Add BoundaryFlows to face flows - call this%addBoundaryFlows_cellRect(cellDefn) + call this%addBoundaryFlows_cellRect(defn) ! -- Set inoexitface flag - cellDefn%inoexitface = 1 + defn%inoexitface = 1 do m = 1, npolyverts + 3 ! kluge note: can be streamlined with above code - if (cellDefn%faceflow(m) < 0d0) cellDefn%inoexitface = 0 + if (defn%faceflow(m) < 0d0) defn%inoexitface = 0 end do ! -- Add up net distributed flow - cellDefn%distflow = this%fmi%SourceFlows(ic) + this%fmi%SinkFlows(ic) + & + defn%distflow = this%fmi%SourceFlows(ic) + this%fmi%SinkFlows(ic) + & this%fmi%StorageFlows(ic) ! -- Set weak sink flag if (this%fmi%SinkFlows(ic) .ne. 0d0) then - cellDefn%iweaksink = 1 + defn%iweaksink = 1 else - cellDefn%iweaksink = 0 + defn%iweaksink = 0 end if end subroutine load_cellDefn_flows !> @brief Add boundary flows to the cell definition faceflow array. !! Assumes cell index and number of vertices are already loaded. - subroutine addBoundaryFlows_cellRect(this, cellDefn) + subroutine addBoundaryFlows_cellRect(this, defn) ! -- dummy class(MethodDisType), intent(inout) :: this - type(CellDefnType), intent(inout) :: cellDefn + type(CellDefnType), intent(inout) :: defn ! -- local integer :: ic, npolyverts integer :: ioffset - ic = cellDefn%icell - npolyverts = cellDefn%npolyverts + ic = defn%icell + npolyverts = defn%npolyverts ioffset = (ic - 1) * 10 - cellDefn%faceflow(1) = cellDefn%faceflow(1) + & + defn%faceflow(1) = defn%faceflow(1) + & this%fmi%BoundaryFlows(ioffset + 1) ! kluge note: should these be additive (seems so)??? - cellDefn%faceflow(2) = cellDefn%faceflow(2) + & + defn%faceflow(2) = defn%faceflow(2) + & this%fmi%BoundaryFlows(ioffset + 2) - cellDefn%faceflow(3) = cellDefn%faceflow(3) + & + defn%faceflow(3) = defn%faceflow(3) + & this%fmi%BoundaryFlows(ioffset + 3) - cellDefn%faceflow(4) = cellDefn%faceflow(4) + & + defn%faceflow(4) = defn%faceflow(4) + & this%fmi%BoundaryFlows(ioffset + 4) - cellDefn%faceflow(5) = cellDefn%faceflow(1) - cellDefn%faceflow(6) = cellDefn%faceflow(6) + & + defn%faceflow(5) = defn%faceflow(1) + defn%faceflow(6) = defn%faceflow(6) + & this%fmi%BoundaryFlows(ioffset + 9) - cellDefn%faceflow(7) = cellDefn%faceflow(7) + & + defn%faceflow(7) = defn%faceflow(7) + & this%fmi%BoundaryFlows(ioffset + 10) end subroutine addBoundaryFlows_cellRect @@ -476,22 +476,22 @@ end subroutine addBoundaryFlows_cellRect !> @brief Load 180-degree vertex indicator array and set flags !! indicating how cell can be represented (kluge: latter needed?). !! Assumes cell index and number of vertices are already loaded. - subroutine load_cellDefn_ispv180(this, cellDefn) + subroutine load_cellDefn_ispv180(this, defn) ! -- dummy class(MethodDisType), intent(inout) :: this - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- allocate ispv180 array - call ExpandArray(cellDefn%ispv180, cellDefn%npolyverts + 1) + call ExpandArray(defn%ispv180, defn%npolyverts + 1) ! -- Load 180-degree indicator. Note that the ispv180 array ! -- does not get reallocated if it is already allocated ! -- to a size greater than or equal to npolyverts+1. - cellDefn%ispv180(1:cellDefn%npolyverts + 1) = .false. + defn%ispv180(1:defn%npolyverts + 1) = .false. ! -- Set flags that indicate how cell can be represented. - cellDefn%canBeCellRect = .true. - cellDefn%canBeCellRectQuad = .false. + defn%canBeCellRect = .true. + defn%canBeCellRectQuad = .false. end subroutine load_cellDefn_ispv180 end module MethodDisModule diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90 index 7e146e43ea8..cab605e2185 100644 --- a/src/Solution/ParticleTracker/MethodDisv.f90 +++ b/src/Solution/ParticleTracker/MethodDisv.f90 @@ -109,21 +109,21 @@ subroutine loadsub_mGDv(this, particle, levelNext, submethod) ! load cell definition ic = particle%iTrackingDomain(levelNext) ! kluge note: is cell number always known coming in? - call this%load_cellDefn(ic, this%cellPoly%cellDefn) + call this%load_cellDefn(ic, this%cellPoly%defn) if (this%fmi%ibdgwfsat0(ic) == 0) then ! kluge note: use cellDefn%sat == DZERO here instead? ! -- Cell is active but dry, so select and initialize pass-to-bottom ! -- cell method and set cell method pointer - call methodCellPassToBot%init(particle, this%cellPoly%cellDefn, & + call methodCellPassToBot%init(particle, this%cellPoly%defn, & this%trackctl) submethod => methodCellPassToBot else ! -- Select and initialize cell method and set cell method pointer - if (this%cellPoly%cellDefn%canBeCellRect) then + if (this%cellPoly%defn%canBeCellRect) then call CellPolyToCellRect(this%cellPoly, cellRect, istatus) call methodCellPollock%init(particle, cellRect, this%trackctl) submethod => methodCellPollock - else if (this%cellPoly%cellDefn%canBeCellRectQuad) then + else if (this%cellPoly%defn%canBeCellRectQuad) then call CellPolyToCellRectQuad(this%cellPoly, cellRectQuad, istatus) call methodCellPollockQuad%init(particle, cellRectQuad, this%trackctl) submethod => methodCellPollockQuad @@ -150,7 +150,7 @@ subroutine pass_mGDv(this, particle) inface = particle%iTrackingDomainBoundary(2) z = particle%z - inbr = this%cellPoly%cellDefn%facenbr(inface) + inbr = this%cellPoly%defn%facenbr(inface) if (inbr .eq. 0) then ! -- Exterior face; no neighbor to map to ! particle%iTrackingDomain(1) = 0 @@ -163,7 +163,7 @@ subroutine pass_mGDv(this, particle) kstp=kstp, reason=3) ! reason=3: termination ! particle%iTrackingDomainBoundary(2) = -1 else - idiag = this%fmi%dis%con%ia(this%cellPoly%cellDefn%icell) + idiag = this%fmi%dis%con%ia(this%cellPoly%defn%icell) ipos = idiag + inbr ic = this%fmi%dis%con%ja(ipos) ! kluge note, todo: use PRT model's DIS instead of fmi's?? particle%iTrackingDomain(2) = ic @@ -177,7 +177,7 @@ subroutine pass_mGDv(this, particle) particle%ilay = ilay end select - call this%mapToNbrCell(this%cellPoly%cellDefn, inface, z) + call this%mapToNbrCell(this%cellPoly%defn, inface, z) particle%iTrackingDomainBoundary(2) = inface particle%iTrackingDomain(3:) = 0 particle%iTrackingDomainBoundary(3:) = 0 @@ -194,33 +194,33 @@ subroutine pass_mGDv(this, particle) end subroutine pass_mGDv !> @brief Map location on cell face to shared cell face of neighbor - subroutine mapToNbrCell(this, cellDefnin, inface, z) + subroutine mapToNbrCell(this, defn, inface, z) ! dummy class(MethodDisvType), intent(inout) :: this - type(CellDefnType), pointer, intent(inout) :: cellDefnin + type(CellDefnType), pointer, intent(inout) :: defn integer, intent(inout) :: inface double precision, intent(inout) :: z ! local integer :: icin, npolyvertsin integer :: ic, npolyverts, inbr, inbrnbr, j, m real(DP) :: zrel, topfrom, botfrom, top, bot, sat - type(CellDefnType), pointer :: cellDefn + type(CellDefnType), pointer :: cd ! -- Map to shared cell face of neighbor - inbr = cellDefnin%facenbr(inface) + inbr = defn%facenbr(inface) if (inbr .eq. 0) then ! kluge note: redundant check ! -- Exterior face; no neighbor to map to inface = -1 ! kluge??? else ! -- Load definition for neighbor cell (neighbor with shared face) - icin = cellDefnin%icell + icin = defn%icell j = this%fmi%dis%con%ia(icin) ic = this%fmi%dis%con%ja(j + inbr) - call create_cellDefn(cellDefn) + call create_cellDefn(cd) ! kluge note: really only need to load facenbr and npolyverts for this - call this%load_cellDefn(ic, cellDefn) ! kluge - npolyvertsin = cellDefnin%npolyverts - npolyverts = cellDefn%npolyverts + call this%load_cellDefn(ic, cd) ! kluge + npolyvertsin = defn%npolyverts + npolyverts = cd%npolyverts if (inface .eq. npolyvertsin + 2) then ! -- Exits through bot, enters through top inface = npolyverts + 3 @@ -232,15 +232,15 @@ subroutine mapToNbrCell(this, cellDefnin, inface, z) j = this%fmi%dis%con%ia(ic) ! kluge note: use shared_edge in DisvGeom to find shared polygon face??? do m = 1, npolyverts + 3 - inbrnbr = cellDefn%facenbr(m) + inbrnbr = cd%facenbr(m) if (this%fmi%dis%con%ja(j + inbrnbr) .eq. icin) then inface = m exit end if end do ! -- Map z between cells - topfrom = cellDefnin%top - botfrom = cellDefnin%bot + topfrom = defn%top + botfrom = defn%bot zrel = (z - botfrom) / (topfrom - botfrom) ! kluge note: use PRT model's DIS instead of fmi's??? top = this%fmi%dis%top(ic) @@ -248,7 +248,7 @@ subroutine mapToNbrCell(this, cellDefnin, inface, z) sat = this%fmi%gwfsat(ic) z = bot + zrel * sat * (top - bot) end if - deallocate (cellDefn) + deallocate (cd) end if end subroutine mapToNbrCell @@ -308,108 +308,108 @@ function get_top(this, iatop) result(top) end function get_top !> @brief Load cell definition from the grid - subroutine load_cellDefn(this, ic, cellDefn) + subroutine load_cellDefn(this, ic, defn) implicit none ! -- dummy class(MethodDisvType), intent(inout) :: this integer, intent(in) :: ic - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- Load basic cell definition components - call this%load_cellDefn_basic(ic, cellDefn) + call this%load_cellDefn_basic(ic, defn) ! -- Load polygon vertices - call this%load_cellDefn_polyverts(cellDefn) + call this%load_cellDefn_polyverts(defn) ! -- Load face neighbors - call this%load_cellDefn_facenbr(cellDefn) + call this%load_cellDefn_facenbr(defn) ! -- Load 180-degree indicator - call this%load_cellDefn_ispv180(cellDefn) + call this%load_cellDefn_ispv180(defn) ! -- Load flows (assumes face neighbors already loaded) - call this%load_cellDefn_flows(cellDefn) + call this%load_cellDefn_flows(defn) end subroutine load_cellDefn !> @brief Loads basic cell definition components from the grid !! These include cell index, vertex count, (ia)top and botm, !! porosity, retfactor, and izone. - subroutine load_cellDefn_basic(this, ic, cellDefn) + subroutine load_cellDefn_basic(this, ic, defn) implicit none ! -- dummy class(MethodDisvType), intent(inout) :: this integer, intent(in) :: ic - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- local ! integer :: iatop,nnbrs integer :: iatop real(DP) :: top, bot, sat ! -- cell index - cellDefn%icell = ic + defn%icell = ic ! -- number of polygon vertices - cellDefn%npolyverts = this%get_npolyverts(ic) + defn%npolyverts = this%get_npolyverts(ic) ! -- iatop, top, and bot iatop = get_iatop(this%fmi%dis%get_ncpl(), & this%fmi%dis%get_nodeuser(ic)) - cellDefn%iatop = iatop + defn%iatop = iatop top = this%fmi%dis%top(ic) bot = this%fmi%dis%bot(ic) sat = this%fmi%gwfsat(ic) top = bot + sat * (top - bot) - cellDefn%top = top - cellDefn%bot = bot - cellDefn%sat = sat + defn%top = top + defn%bot = bot + defn%sat = sat ! -- Load porosity, retfactor, and izone - cellDefn%porosity = this%porosity(ic) - cellDefn%retfactor = this%retfactor(ic) - cellDefn%izone = this%izone(ic) + defn%porosity = this%porosity(ic) + defn%retfactor = this%retfactor(ic) + defn%izone = this%izone(ic) end subroutine load_cellDefn_basic !> @brief Loads polygon vertices to cell definition from the grid !! Assumes cell index is already loaded into cell definition. - subroutine load_cellDefn_polyverts(this, cellDefn) + subroutine load_cellDefn_polyverts(this, defn) ! -- dummy class(MethodDisvType), intent(inout) :: this - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- Load polygon vertices call this%fmi%dis%get_polyverts( & - cellDefn%icell, & - cellDefn%polyvert, & + defn%icell, & + defn%polyvert, & closed=.true.) end subroutine load_cellDefn_polyverts !> @brief Loads face neighbors to cell definition from the grid !! Assumes cell index and number of vertices are already loaded. - subroutine load_cellDefn_facenbr(this, cellDefn) + subroutine load_cellDefn_facenbr(this, defn) use GwfDisvModule implicit none ! -- dummy class(MethodDisvType), intent(inout) :: this - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- local integer :: ic, npolyverts integer :: ic1, ic2, icu1, icu2, j1, j2, k1, k2, iloc, ipos integer :: istart1, istart2, istop1, istop2, iedgeface integer :: ncpl - ic = cellDefn%icell - npolyverts = cellDefn%npolyverts + ic = defn%icell + npolyverts = defn%npolyverts ! -- Load face neighbors. Note that the facenbr array ! -- does not get reallocated if it is already allocated ! -- to a size greater than or equal to npolyverts+3. - call allocate_as_needed(cellDefn%facenbr, npolyverts + 3) + call allocate_as_needed(defn%facenbr, npolyverts + 3) select type (dis => this%fmi%dis) ! kluge type guard type is (GwfDisvType) ! kluge - cellDefn%facenbr = 0 + defn%facenbr = 0 ic1 = ic icu1 = dis%get_nodeuser(ic1) ncpl = dis%get_ncpl() @@ -429,21 +429,21 @@ subroutine load_cellDefn_facenbr(this, cellDefn) iedgeface) if (iedgeface /= 0) then ! -- Edge (polygon) face neighbor - cellDefn%facenbr(iedgeface) = int(iloc, 1) + defn%facenbr(iedgeface) = int(iloc, 1) else if (k2 > k1) then ! -- Bottom face neighbor - cellDefn%facenbr(npolyverts + 2) = int(iloc, 1) + defn%facenbr(npolyverts + 2) = int(iloc, 1) else if (k2 < k1) then ! -- Top face neighbor - cellDefn%facenbr(npolyverts + 3) = int(iloc, 1) + defn%facenbr(npolyverts + 3) = int(iloc, 1) else print *, "Error -- k2 should be <> k1, since no shared edge face" end if end if end do ! -- List of edge (polygon) faces wraps around - cellDefn%facenbr(npolyverts + 1) = cellDefn%facenbr(1) + defn%facenbr(npolyverts + 1) = defn%facenbr(1) end select end subroutine load_cellDefn_facenbr @@ -489,98 +489,98 @@ end subroutine shared_edgeface !> @brief Load flows into the cell definition. !! These include face flows and net distributed flows. !! Assumes cell index and number of vertices are already loaded. - subroutine load_cellDefn_flows(this, cellDefn) + subroutine load_cellDefn_flows(this, defn) ! -- dummy class(MethodDisvType), intent(inout) :: this - type(CellDefnType), intent(inout) :: cellDefn + type(CellDefnType), intent(inout) :: defn ! -- local integer :: ic, npolyverts, m, n - ic = cellDefn%icell - npolyverts = cellDefn%npolyverts + ic = defn%icell + npolyverts = defn%npolyverts ! -- Load face flows. Note that the faceflow array ! -- does not get reallocated if it is already allocated ! -- to a size greater than or equal to npolyverts+3. - call allocate_as_needed(cellDefn%faceflow, npolyverts + 3) - cellDefn%faceflow = 0d0 + call allocate_as_needed(defn%faceflow, npolyverts + 3) + defn%faceflow = 0d0 ! -- As with polygon nbrs, polygon face flows wrap around for ! -- convenience at position npolyverts+1, and bot and top flows ! -- are tacked on the end of the list do m = 1, npolyverts + 3 - n = cellDefn%facenbr(m) + n = defn%facenbr(m) if (n > 0) & - cellDefn%faceflow(m) = this%fmi%gwfflowja(this%fmi%dis%con%ia(ic) + n) + defn%faceflow(m) = this%fmi%gwfflowja(this%fmi%dis%con%ia(ic) + n) end do - call this%addBoundaryFlows_cellPoly(cellDefn) + call this%addBoundaryFlows_cellPoly(defn) ! -- Set inoexitface flag - cellDefn%inoexitface = 1 + defn%inoexitface = 1 do m = 1, npolyverts + 3 ! kluge note: can be streamlined with above code - if (cellDefn%faceflow(m) < 0d0) cellDefn%inoexitface = 0 + if (defn%faceflow(m) < 0d0) defn%inoexitface = 0 end do ! -- Add up net distributed flow - cellDefn%distflow = this%fmi%SourceFlows(ic) + this%fmi%SinkFlows(ic) + & + defn%distflow = this%fmi%SourceFlows(ic) + this%fmi%SinkFlows(ic) + & this%fmi%StorageFlows(ic) ! -- Set weak sink flag if (this%fmi%SinkFlows(ic) .ne. 0d0) then - cellDefn%iweaksink = 1 + defn%iweaksink = 1 else - cellDefn%iweaksink = 0 + defn%iweaksink = 0 end if end subroutine load_cellDefn_flows !> @brief Load boundary flows from the grid into a rectangular cell. !! Assumes cell index and number of vertices are already loaded. - subroutine addBoundaryFlows_cellRect(this, cellDefn) + subroutine addBoundaryFlows_cellRect(this, defn) ! -- dummy class(MethodDisvType), intent(inout) :: this - type(CellDefnType), intent(inout) :: cellDefn + type(CellDefnType), intent(inout) :: defn ! -- local integer :: ic, npolyverts integer :: ioffset - ic = cellDefn%icell - npolyverts = cellDefn%npolyverts + ic = defn%icell + npolyverts = defn%npolyverts ! kluge note - assignment of BoundaryFlows to faceflow below assumes vertex 1 ! is at upper left of rectangular cell, and BoundaryFlows use old iface order ! ioffset = (ic - 1)*6 ioffset = (ic - 1) * 10 ! kluge note: should these be additive (seems so)??? - cellDefn%faceflow(1) = cellDefn%faceflow(1) + & + defn%faceflow(1) = defn%faceflow(1) + & this%fmi%BoundaryFlows(ioffset + 4) - cellDefn%faceflow(2) = cellDefn%faceflow(2) + & + defn%faceflow(2) = defn%faceflow(2) + & this%fmi%BoundaryFlows(ioffset + 2) - cellDefn%faceflow(3) = cellDefn%faceflow(3) + & + defn%faceflow(3) = defn%faceflow(3) + & this%fmi%BoundaryFlows(ioffset + 3) - cellDefn%faceflow(4) = cellDefn%faceflow(4) + & + defn%faceflow(4) = defn%faceflow(4) + & this%fmi%BoundaryFlows(ioffset + 1) - cellDefn%faceflow(5) = cellDefn%faceflow(1) - cellDefn%faceflow(6) = cellDefn%faceflow(6) + & + defn%faceflow(5) = defn%faceflow(1) + defn%faceflow(6) = defn%faceflow(6) + & this%fmi%BoundaryFlows(ioffset + 9) - cellDefn%faceflow(7) = cellDefn%faceflow(7) + & + defn%faceflow(7) = defn%faceflow(7) + & this%fmi%BoundaryFlows(ioffset + 10) end subroutine addBoundaryFlows_cellRect !> @brief Load boundary flows from the grid into rectangular quadcell. !! Assumes cell index and number of vertices are already loaded. - subroutine addBoundaryFlows_cellRectQuad(this, cellDefn) + subroutine addBoundaryFlows_cellRectQuad(this, defn) ! -- dummy class(MethodDisvType), intent(inout) :: this - type(CellDefnType), intent(inout) :: cellDefn + type(CellDefnType), intent(inout) :: defn ! -- local integer :: ic, npolyverts, m, n, nn integer :: ioffset, nbf, m1, m2, mdiff double precision :: qbf integer :: irectvert(5) ! kluge - ic = cellDefn%icell - npolyverts = cellDefn%npolyverts + ic = defn%icell + npolyverts = defn%npolyverts ! kluge note - assignment of BoundaryFlows to faceflow below assumes vertex 1 ! is at upper left of rectangular cell, and BoundaryFlows use old iface order @@ -598,7 +598,7 @@ subroutine addBoundaryFlows_cellRectQuad(this, cellDefn) qbf = this%fmi%BoundaryFlows(ioffset + nbf) nn = 0 ! kluge ... do m = 1, npolyverts - if (.not. cellDefn%ispv180(m)) then + if (.not. defn%ispv180(m)) then nn = nn + 1 irectvert(nn) = m end if @@ -610,54 +610,54 @@ subroutine addBoundaryFlows_cellRectQuad(this, cellDefn) mdiff = m2 - m1 if (mdiff .eq. 1) then ! -- Assign BoundaryFlow to corresponding polygon face - cellDefn%faceflow(m1) = cellDefn%faceflow(m1) + qbf + defn%faceflow(m1) = defn%faceflow(m1) + qbf else ! -- Split BoundaryFlow between two faces on quad-refined edge qbf = 5d-1 * qbf - cellDefn%faceflow(m1) = cellDefn%faceflow(m1) + qbf - cellDefn%faceflow(m1 + 1) = cellDefn%faceflow(m1 + 1) + qbf + defn%faceflow(m1) = defn%faceflow(m1) + qbf + defn%faceflow(m1 + 1) = defn%faceflow(m1 + 1) + qbf end if end do ! -- Wrap around to 1 in position npolyverts+1 m = npolyverts + 1 - cellDefn%faceflow(m) = cellDefn%faceflow(1) + defn%faceflow(m) = defn%faceflow(1) ! -- Bottom in position npolyverts+2 m = m + 1 - cellDefn%faceflow(m) = cellDefn%faceflow(m) + & + defn%faceflow(m) = defn%faceflow(m) + & this%fmi%BoundaryFlows(ioffset + 9) ! -- Top in position npolyverts+3 m = m + 1 - cellDefn%faceflow(m) = cellDefn%faceflow(m) + & + defn%faceflow(m) = defn%faceflow(m) + & this%fmi%BoundaryFlows(ioffset + 10) end subroutine addBoundaryFlows_cellRectQuad !> @brief Load boundary flows from the grid into a polygonal cell. !! Assumes cell index and number of vertices are already loaded. - subroutine addBoundaryFlows_cellPoly(this, cellDefn) + subroutine addBoundaryFlows_cellPoly(this, defn) ! -- dummy class(MethodDisvType), intent(inout) :: this - type(CellDefnType), intent(inout) :: cellDefn + type(CellDefnType), intent(inout) :: defn ! -- local integer :: ic, npolyverts, ioffset, iv - ic = cellDefn%icell - npolyverts = cellDefn%npolyverts + ic = defn%icell + npolyverts = defn%npolyverts ! kluge note: hardwired for max 8 polygon faces plus top and bottom for now ioffset = (ic - 1) * 10 do iv = 1, npolyverts ! kluge note: should these be additive (seems so)??? - cellDefn%faceflow(iv) = & - cellDefn%faceflow(iv) + & + defn%faceflow(iv) = & + defn%faceflow(iv) + & this%fmi%BoundaryFlows(ioffset + iv) end do - cellDefn%faceflow(npolyverts + 1) = cellDefn%faceflow(1) - cellDefn%faceflow(npolyverts + 2) = & - cellDefn%faceflow(npolyverts + 2) + & + defn%faceflow(npolyverts + 1) = defn%faceflow(1) + defn%faceflow(npolyverts + 2) = & + defn%faceflow(npolyverts + 2) + & this%fmi%BoundaryFlows(ioffset + 9) - cellDefn%faceflow(npolyverts + 3) = & - cellDefn%faceflow(npolyverts + 3) + & + defn%faceflow(npolyverts + 3) = & + defn%faceflow(npolyverts + 3) + & this%fmi%BoundaryFlows(ioffset + 10) end subroutine addBoundaryFlows_cellPoly @@ -665,10 +665,10 @@ end subroutine addBoundaryFlows_cellPoly !> @brief Load 180-degree vertex indicator array and set flags !! indicating how cell can be represented (kluge: latter needed?). !! Assumes cell index and number of vertices are already loaded. - subroutine load_cellDefn_ispv180(this, cellDefn) ! kluge note: rename??? + subroutine load_cellDefn_ispv180(this, defn) ! kluge note: rename??? ! -- dummy class(MethodDisvType), intent(inout) :: this - type(CellDefnType), pointer, intent(inout) :: cellDefn + type(CellDefnType), pointer, intent(inout) :: defn ! -- local integer :: npolyverts, m, m0, m1, m2 integer :: ic @@ -678,17 +678,17 @@ subroutine load_cellDefn_ispv180(this, cellDefn) ! kluge note: rename??? s0mag, s2x, s2y, s2mag, sinang, dotprod logical last180 - ic = cellDefn%icell - npolyverts = cellDefn%npolyverts + ic = defn%icell + npolyverts = defn%npolyverts ! -- Load 180-degree indicator. Note that the ispv180 array ! -- does not get reallocated if it is already allocated ! -- to a size greater than or equal to npolyverts+1. ! -- Also, set flags that indicate how cell can be represented. - call allocate_as_needed(cellDefn%ispv180, npolyverts + 1) - cellDefn%ispv180(1:npolyverts + 1) = .false. - cellDefn%canBeCellRect = .false. - cellDefn%canBeCellRectQuad = .false. + call allocate_as_needed(defn%ispv180, npolyverts + 1) + defn%ispv180(1:npolyverts + 1) = .false. + defn%canBeCellRect = .false. + defn%canBeCellRectQuad = .false. epsang = 1d-3 ! kluge hardwire, and using one value for all angles epslen = 1d-3 ! kluge hardwire num90 = 0 @@ -709,12 +709,12 @@ subroutine load_cellDefn_ispv180(this, cellDefn) ! kluge note: rename??? m0 = m1 - 1 m2 = m1 + 1 end if - x0 = cellDefn%polyvert(1, m0) - y0 = cellDefn%polyvert(2, m0) - x1 = cellDefn%polyvert(1, m1) - y1 = cellDefn%polyvert(2, m1) - x2 = cellDefn%polyvert(1, m2) - y2 = cellDefn%polyvert(2, m2) + x0 = defn%polyvert(1, m0) + y0 = defn%polyvert(2, m0) + x1 = defn%polyvert(1, m1) + y1 = defn%polyvert(2, m1) + x2 = defn%polyvert(1, m2) + y2 = defn%polyvert(2, m2) s0x = x0 - x1 s0y = y0 - y1 s0mag = dsqrt(s0x * s0x + s0y * s0y) @@ -748,7 +748,7 @@ subroutine load_cellDefn_ispv180(this, cellDefn) ! kluge note: rename??? ! one place (procedure) to avoid potential disparities between multiple checks num180 = num180 + 1 last180 = .true. - cellDefn%ispv180(m) = .true. + defn%ispv180(m) = .true. end if else if (sinang .gt. 0d0) then numacute = numacute + 1 @@ -770,15 +770,15 @@ subroutine load_cellDefn_ispv180(this, cellDefn) ! kluge note: rename??? stop end if ! -- List of 180-degree indicators wraps around for convenience - cellDefn%ispv180(npolyverts + 1) = cellDefn%ispv180(1) + defn%ispv180(npolyverts + 1) = defn%ispv180(1) ! if (num90 .eq. 4) then if (num180 .eq. 0) then ! -- Can become CellRect - cellDefn%canBeCellRect = .true. + defn%canBeCellRect = .true. else ! -- Can become CellRectQuad - cellDefn%canBeCellRectQuad = .true. + defn%canBeCellRectQuad = .true. end if end if