disvgeom Module Reference

Data Types
type	disvgeomtype

Functions/Subroutines
subroutine	init (this, nlay, ncpl, nodes, top_grid, bot_grid, iavert, javert, vertex_grid, cellxy_grid, nodereduced, nodeuser)
	Initialize. More...
subroutine	set_kj (this, k, j)
	Set node IDs. More...
subroutine	set_nodered (this, nodered)
	Set reduced node number. More...
subroutine	cell_setup (this)
	Set top and bottom elevations of grid cell. More...
subroutine	cprops (this, cell2, hwva, cl1, cl2, ax, ihc)
subroutine	edge_normal (this, cell2, xcomp, ycomp)
	Return the x and y components of an outward normal facing vector. More...
subroutine	connection_vector (this, cell2, nozee, satn, satm, xcomp, ycomp, zcomp, conlen)
	Return the x y and z components of a unit vector that points from from the center of this to the center of cell2, and the straight-line connection length. More...
logical function	shares_edge (this, cell2)
	Return true if this shares a horizontal edge with cell2. More...
subroutine, public	shared_edge (ivlist1, ivlist2, ivert1, ivert2)
	Find two common vertices shared by cell1 and cell2. More...
real(dp) function	get_area (this)
	Calculate and return the area of the cell. More...
real(dp) function	anglex (x1, y1, x2, y2)
	Calculate the angle that the x-axis makes with a line that is normal to the two points. More...
real(dp) function	distance (x1, y1, x2, y2)
	Calculate distance between two points. More...
real(dp) function	distance_normal (x0, y0, x1, y1, x2, y2)
	Calculate normal distance from point (x0, y0) to line defined by two points, (x1, y1), (x2, y2). More...
subroutine	line_unit_normal (x0, y0, x1, y1, xcomp, ycomp)
	Calculate the normal vector components (xcomp and ycomp) for a line defined by two points, (x0, y0), (x1, y1). More...
subroutine, public	line_unit_vector (x0, y0, z0, x1, y1, z1, xcomp, ycomp, zcomp, vmag)
	Calculate the vector components (xcomp, ycomp, and zcomp) for a line defined by two points, (x0, y0, z0), (x1, y1, z1). More...

Function/Subroutine Documentation

anglex()

real(dp) function disvgeom::anglex	(	real(dp), intent(in)	x1,
		real(dp), intent(in)	y1,
		real(dp), intent(in)	x2,
		real(dp), intent(in)	y2
	)

This assumes that vertices are numbered clockwise so that the angle is for the normal outward of cell n.

Definition at line 394 of file DisvGeom.f90.

    ! -- modules
    use constantsmodule, only: dzero, dtwo, dpi
    ! -- dummy
    real(DP), intent(in) :: x1
    real(DP), intent(in) :: x2
    real(DP), intent(in) :: y1
    real(DP), intent(in) :: y2
    ! -- return
    real(DP) :: ax
    ! -- local
    real(DP) :: dx
    real(DP) :: dy
    !
    dx = x2 - x1
    dy = y2 - y1
    ax = atan2(dx, -dy)
    if (ax < dzero) ax = dtwo * dpi + ax

Here is the caller graph for this function:

cell_setup()

subroutine disvgeom::cell_setup ( class(disvgeomtype)  this )

private

Definition at line 128 of file DisvGeom.f90.

    ! -- dummy
    class(DisvGeomType) :: this
    !
    this%top = this%top_grid(this%nodered)
    this%bot = this%bot_grid(this%nodered)

connection_vector()

subroutine disvgeom::connection_vector	(	class(disvgeomtype)	this,
		type(disvgeomtype)	cell2,
		logical, intent(in)	nozee,
		real(dp), intent(in)	satn,
		real(dp), intent(in)	satm,
		real(dp), intent(out)	xcomp,
		real(dp), intent(out)	ycomp,
		real(dp), intent(out)	zcomp,
		real(dp), intent(out)	conlen
	)

Definition at line 238 of file DisvGeom.f90.

    ! -- module
    use constantsmodule, only: dzero, dhalf, done
    ! -- dummy
    class(DisvGeomType) :: this
    type(DisvGeomType) :: cell2
    logical, intent(in) :: nozee
    real(DP), intent(in) :: satn
    real(DP), intent(in) :: satm
    real(DP), intent(out) :: xcomp
    real(DP), intent(out) :: ycomp
    real(DP), intent(out) :: zcomp
    real(DP), intent(out) :: conlen
    ! -- local
    real(DP) :: x1, y1, z1, x2, y2, z2
    !
    x1 = this%cellxy_grid(1, this%j)
    y1 = this%cellxy_grid(2, this%j)
    x2 = this%cellxy_grid(1, cell2%j)
    y2 = this%cellxy_grid(2, cell2%j)
    if (nozee) then
      z1 = dzero
      z2 = dzero
    else
      z1 = this%bot + dhalf * satn * (this%top - this%bot)
      z2 = cell2%bot + dhalf * satm * (cell2%top - cell2%bot)
    end if
    !
    call line_unit_vector(x1, y1, z1, x2, y2, z2, xcomp, ycomp, zcomp, &
                          conlen)

Here is the call graph for this function:

cprops()

subroutine disvgeom::cprops ( class(disvgeomtype)  this, type(disvgeomtype)  cell2, real(dp), intent(out)  hwva, real(dp), intent(out)  cl1, real(dp), intent(out)  cl2, real(dp), intent(out)  ax, integer(i4b), intent(out)  ihc  )

private

Definition at line 136 of file DisvGeom.f90.

    ! -- module
    use constantsmodule, only: dzero, dhalf, done
    ! -- dummy
    class(DisvGeomType) :: this
    type(DisvGeomType) :: cell2
    real(DP), intent(out) :: hwva
    real(DP), intent(out) :: cl1
    real(DP), intent(out) :: cl2
    real(DP), intent(out) :: ax
    integer(I4B), intent(out) :: ihc
    ! -- local
    integer(I4B) :: ivert1, ivert2
    integer(I4B) :: istart1, istart2, istop1, istop2
    real(DP) :: x0, y0, x1, y1, x2, y2
    !
    if (this%j == cell2%j) then
      !
      ! -- Cells share same j index, so must be a vertical connection
      ihc = 0
      hwva = this%get_area()
      cl1 = dhalf * (this%top - this%bot)
      cl2 = dhalf * (cell2%top - cell2%bot)
      ax = dzero
    else
      !
      ! -- Must be horizontal connection
      ihc = 1
      istart1 = this%iavert(this%j)
      istop1 = this%iavert(this%j + 1) - 1
      istart2 = cell2%iavert(cell2%j)
      istop2 = this%iavert(cell2%j + 1) - 1
      call shared_edge(this%javert(istart1:istop1), &
                       this%javert(istart2:istop2), &
                       ivert1, ivert2)
      if (ivert1 == 0 .or. ivert2 == 0) then
        !
        ! -- Cells do not share an edge
        hwva = dzero
        cl1 = done
        cl2 = done
      else
        x1 = this%vertex_grid(1, ivert1)
        y1 = this%vertex_grid(2, ivert1)
        x2 = this%vertex_grid(1, ivert2)
        y2 = this%vertex_grid(2, ivert2)
        hwva = distance(x1, y1, x2, y2)
        !
        ! -- cl1
        x0 = this%cellxy_grid(1, this%j)
        y0 = this%cellxy_grid(2, this%j)
        cl1 = distance_normal(x0, y0, x1, y1, x2, y2)
        !
        ! -- cl2
        x0 = this%cellxy_grid(1, cell2%j)
        y0 = this%cellxy_grid(2, cell2%j)
        cl2 = distance_normal(x0, y0, x1, y1, x2, y2)
        !
        ! -- anglex
        x1 = this%vertex_grid(1, ivert1)
        y1 = this%vertex_grid(2, ivert1)
        x2 = this%vertex_grid(1, ivert2)
        y2 = this%vertex_grid(2, ivert2)
        ax = anglex(x1, y1, x2, y2)
      end if
    end if

Here is the call graph for this function:

distance()

real(dp) function disvgeom::distance	(	real(dp), intent(in)	x1,
		real(dp), intent(in)	y1,
		real(dp), intent(in)	x2,
		real(dp), intent(in)	y2
	)

Definition at line 416 of file DisvGeom.f90.

    ! -- dummy
    real(DP), intent(in) :: x1
    real(DP), intent(in) :: x2
    real(DP), intent(in) :: y1
    real(DP), intent(in) :: y2
    ! -- return
    real(DP) :: d
    !
    d = (x1 - x2)**2 + (y1 - y2)**2
    d = sqrt(d)

Here is the caller graph for this function:

distance_normal()

real(dp) function disvgeom::distance_normal ( real(dp), intent(in)  x0, real(dp), intent(in)  y0, real(dp), intent(in)  x1, real(dp), intent(in)  y1, real(dp), intent(in)  x2, real(dp), intent(in)  y2  )

private

Definition at line 432 of file DisvGeom.f90.

    ! -- dummy
    real(DP), intent(in) :: x0
    real(DP), intent(in) :: y0
    real(DP), intent(in) :: x1
    real(DP), intent(in) :: y1
    real(DP), intent(in) :: x2
    real(DP), intent(in) :: y2
    ! -- return
    real(DP) :: d
    !
    d = abs((x2 - x1) * (y1 - y0) - (x1 - x0) * (y2 - y1))
    d = d / distance(x1, y1, x2, y2)

Here is the call graph for this function:

Here is the caller graph for this function:

edge_normal()

subroutine disvgeom::edge_normal	(	class(disvgeomtype)	this,
		type(disvgeomtype)	cell2,
		real(dp), intent(out)	xcomp,
		real(dp), intent(out)	ycomp
	)

Definition at line 206 of file DisvGeom.f90.

    ! -- module
    use constantsmodule, only: dzero, dhalf, done
    ! -- dummy
    class(DisvGeomType) :: this
    type(DisvGeomType) :: cell2
    real(DP), intent(out) :: xcomp
    real(DP), intent(out) :: ycomp
    ! -- local
    integer(I4B) :: ivert1, ivert2
    integer(I4B) :: istart1, istart2, istop1, istop2
    real(DP) :: x1, y1, x2, y2
    !
    istart1 = this%iavert(this%j)
    istop1 = this%iavert(this%j + 1) - 1
    istart2 = cell2%iavert(cell2%j)
    istop2 = this%iavert(cell2%j + 1) - 1
    call shared_edge(this%javert(istart1:istop1), &
                     this%javert(istart2:istop2), &
                     ivert1, ivert2)
    x1 = this%vertex_grid(1, ivert1)
    y1 = this%vertex_grid(2, ivert1)
    x2 = this%vertex_grid(1, ivert2)
    y2 = this%vertex_grid(2, ivert2)
    !
    call line_unit_normal(x1, y1, x2, y2, xcomp, ycomp)

Here is the call graph for this function:

get_area()

real(dp) function disvgeom::get_area ( class(disvgeomtype)  this )

private

a = 1/2 *[(x1*y2 + x2*y3 + x3*y4 + ... + xn*y1) - (x2*y1 + x3*y2 + x4*y3 + ... + x1*yn)]

Definition at line 339 of file DisvGeom.f90.

    ! -- module
    use constantsmodule, only: dzero, dhalf
    ! -- dummy
    class(DisvGeomType) :: this
    ! -- return
    real(DP) :: area
    ! -- local
    integer(I4B) :: ivert
    integer(I4B) :: nvert
    integer(I4B) :: icount
    integer(I4B) :: iv1
    real(DP) :: x
    real(DP) :: y
    real(DP) :: x1
    real(DP) :: y1
    !
    area = dzero
    nvert = this%iavert(this%j + 1) - this%iavert(this%j)
    icount = 1
    iv1 = this%javert(this%iavert(this%j))
    x1 = this%vertex_grid(1, iv1)
    y1 = this%vertex_grid(2, iv1)
    do ivert = this%iavert(this%j), this%iavert(this%j + 1) - 1
      x = this%vertex_grid(1, this%javert(ivert))
      if (icount < nvert) then
        y = this%vertex_grid(2, this%javert(ivert + 1))
      else
        y = this%vertex_grid(2, this%javert(this%iavert(this%j)))
      end if
      area = area + (x - x1) * (y - y1)
      icount = icount + 1
    end do
    !
    icount = 1
    do ivert = this%iavert(this%j), this%iavert(this%j + 1) - 1
      y = this%vertex_grid(2, this%javert(ivert))
      if (icount < nvert) then
        x = this%vertex_grid(1, this%javert(ivert + 1))
      else
        x = this%vertex_grid(1, this%javert(this%iavert(this%j)))
      end if
      area = area - (x - x1) * (y - y1)
      icount = icount + 1
    end do
    !
    area = abs(area) * dhalf

init()

subroutine disvgeom::init ( class(disvgeomtype)  this, integer(i4b), intent(in)  nlay, integer(i4b), intent(in)  ncpl, integer(i4b), intent(in)  nodes, real(dp), dimension(nodes), target  top_grid, real(dp), dimension(nodes), target  bot_grid, integer(i4b), dimension(:), target  iavert, integer(i4b), dimension(:), target  javert, real(dp), dimension(:, :), target  vertex_grid, real(dp), dimension(:, :), target  cellxy_grid, integer(i4b), dimension(ncpl, nlay), target  nodereduced, integer(i4b), dimension(nodes), target  nodeuser  )

private

Definition at line 50 of file DisvGeom.f90.

    ! -- dummy
    class(DisvGeomType) :: this
    integer(I4B), intent(in) :: nlay
    integer(I4B), intent(in) :: ncpl
    integer(I4B), intent(in) :: nodes
    real(DP), dimension(nodes), target :: top_grid
    real(DP), dimension(nodes), target :: bot_grid
    integer(I4B), dimension(:), target :: iavert
    integer(I4B), dimension(:), target :: javert
    real(DP), dimension(:, :), target :: vertex_grid
    real(DP), dimension(:, :), target :: cellxy_grid
    integer(I4B), dimension(ncpl, nlay), target :: nodereduced
    integer(I4B), dimension(nodes), target :: nodeuser
    ! -- local
    integer(I4B) :: nodesuser
    !
    this%nlay = nlay
    this%ncpl = ncpl
    this%nodes = nodes
    this%top_grid => top_grid
    this%bot_grid => bot_grid
    this%iavert => iavert
    this%javert => javert
    this%vertex_grid => vertex_grid
    this%cellxy_grid => cellxy_grid
    this%nodereduced => nodereduced
    this%nodeuser => nodeuser
    nodesuser = ncpl * nlay
    !
    if (nodes < nodesuser) then
      this%reduced = .true.
    else
      this%reduced = .false.
    end if

line_unit_normal()

subroutine disvgeom::line_unit_normal ( real(dp), intent(in)  x0, real(dp), intent(in)  y0, real(dp), intent(in)  x1, real(dp), intent(in)  y1, real(dp), intent(out)  xcomp, real(dp), intent(out)  ycomp  )

private

Definition at line 450 of file DisvGeom.f90.

    ! -- dummy
    real(DP), intent(in) :: x0
    real(DP), intent(in) :: y0
    real(DP), intent(in) :: x1
    real(DP), intent(in) :: y1
    real(DP), intent(out) :: xcomp
    real(DP), intent(out) :: ycomp
    ! -- local
    real(DP) :: dx, dy, vmag
    !
    dx = x1 - x0
    dy = y1 - y0
    vmag = sqrt(dx**2 + dy**2)
    xcomp = -dy / vmag
    ycomp = dx / vmag

Here is the caller graph for this function:

line_unit_vector()

subroutine, public disvgeom::line_unit_vector	(	real(dp), intent(in)	x0,
		real(dp), intent(in)	y0,
		real(dp), intent(in)	z0,
		real(dp), intent(in)	x1,
		real(dp), intent(in)	y1,
		real(dp), intent(in)	z1,
		real(dp), intent(out)	xcomp,
		real(dp), intent(out)	ycomp,
		real(dp), intent(out)	zcomp,
		real(dp)	vmag
	)

Also return the magnitude of the original vector, vmag.

Definition at line 473 of file DisvGeom.f90.

    ! -- dummy
    real(DP), intent(in) :: x0
    real(DP), intent(in) :: y0
    real(DP), intent(in) :: z0
    real(DP), intent(in) :: x1
    real(DP), intent(in) :: y1
    real(DP), intent(in) :: z1
    real(DP), intent(out) :: xcomp
    real(DP), intent(out) :: ycomp
    real(DP), intent(out) :: zcomp
    real(DP) :: vmag
    ! -- local
    real(DP) :: dx, dy, dz
    !
    dx = x1 - x0
    dy = y1 - y0
    dz = z1 - z0
    vmag = sqrt(dx**2 + dy**2 + dz**2)
    xcomp = dx / vmag
    ycomp = dy / vmag
    zcomp = dz / vmag

Here is the caller graph for this function:

set_kj()

subroutine disvgeom::set_kj ( class(disvgeomtype)  this, integer(i4b), intent(in)  k, integer(i4b), intent(in)  j  )

private

Definition at line 90 of file DisvGeom.f90.

    ! -- dummy
    class(DisvGeomType) :: this
    integer(I4B), intent(in) :: k
    integer(I4B), intent(in) :: j
    !
    this%k = k
    this%j = j
    this%nodeusr = get_node(k, 1, j, this%nlay, 1, this%ncpl)
    if (this%reduced) then
      this%nodered = this%nodereduced(k, j)
    else
      this%nodered = this%nodeusr
    end if
    call this%cell_setup()

Here is the call graph for this function:

Here is the caller graph for this function:

set_nodered()

subroutine disvgeom::set_nodered ( class(disvgeomtype)  this, integer(i4b), intent(in)  nodered  )

private

Definition at line 109 of file DisvGeom.f90.

    ! -- dummy
    class(DisvGeomType) :: this
    integer(I4B), intent(in) :: nodered
    !
    this%nodered = nodered
    !
    if (this%reduced) then
      this%nodeusr = this%nodeuser(nodered)
    else
      this%nodeusr = nodered
    end if
    !
    call get_jk(this%nodeusr, this%ncpl, this%nlay, this%j, this%k)
    call this%cell_setup()

Here is the call graph for this function:

Here is the caller graph for this function:

shared_edge()

subroutine, public disvgeom::shared_edge	(	integer(i4b), dimension(:)	ivlist1,
		integer(i4b), dimension(:)	ivlist2,
		integer(i4b), intent(out)	ivert1,
		integer(i4b), intent(out)	ivert2
	)

Return 0 if there are no shared edges. Proceed forward through ivlist1 and backward through ivlist2 as a clockwise face in cell1 must correspond to a counter clockwise face in cell2.

Definition at line 302 of file DisvGeom.f90.

    ! -- dummy
    integer(I4B), dimension(:) :: ivlist1
    integer(I4B), dimension(:) :: ivlist2
    integer(I4B), intent(out) :: ivert1
    integer(I4B), intent(out) :: ivert2
    ! -- local
    integer(I4B) :: nv1
    integer(I4B) :: nv2
    integer(I4B) :: il1
    integer(I4B) :: il2
    logical :: found
    !
    found = .false.
    nv1 = size(ivlist1)
    nv2 = size(ivlist2)
    ivert1 = 0
    ivert2 = 0
    outerloop: do il1 = 1, nv1 - 1
      do il2 = nv2, 2, -1
        if (ivlist1(il1) == ivlist2(il2) .and. &
            ivlist1(il1 + 1) == ivlist2(il2 - 1)) then
          found = .true.
          ivert1 = ivlist1(il1)
          ivert2 = ivlist1(il1 + 1)
          exit outerloop
        end if
      end do
      if (found) exit
    end do outerloop

Here is the caller graph for this function:

shares_edge()

logical function disvgeom::shares_edge	(	class(disvgeomtype)	this,
		type(disvgeomtype)	cell2
	)

Definition at line 273 of file DisvGeom.f90.

    ! -- dummy
    class(DisvGeomType) :: this
    type(DisvGeomType) :: cell2
    ! -- return
    logical l
    ! -- local
    integer(I4B) :: istart1, istop1, istart2, istop2
    integer(I4B) :: ivert1, ivert2
    !
    istart1 = this%iavert(this%j)
    istop1 = this%iavert(this%j + 1) - 1
    istart2 = cell2%iavert(cell2%j)
    istop2 = this%iavert(cell2%j + 1) - 1
    call shared_edge(this%javert(istart1:istop1), &
                     this%javert(istart2:istop2), &
                     ivert1, ivert2)
    l = .true.
    if (ivert1 == 0 .or. ivert2 == 0) then
      l = .false.
    end if

Here is the call graph for this function: