Disu.f90

Go to the documentation of this file.

module disumodule
 
  use arrayreadersmodule, only: readarray
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: linelength, lenbigline, lenmempath, lenvarname, &
                             dzero, done, dhalf, dem6, dpio180
  use connectionsmodule, only: iac_to_ia
  use inputoutputmodule, only: urword, ulasav, ulaprufw, ubdsv1, ubdsv06, &
                               getunit, openfile
  use simmodule, only: count_errors, store_error, store_error_unit, &
                       store_error_filename, store_warning
  use simvariablesmodule, only: errmsg, idm_context
  use messagemodule, only: write_message_counter
  use basedismodule, only: disbasetype
  use memorymanagermodule, only: mem_allocate, mem_deallocate, &
                                 mem_reallocate, mem_setptr
  use memorymanagerextmodule, only: mem_set_value, memorystore_remove, &
                                    memorystore_release
  use tdismodule, only: kstp, kper, pertim, totim, delt
  use disvgeom, only: line_unit_vector
 
  implicit none
 
  private
  public :: disutype
  public :: disu_cr
  public :: castasdisutype
 
  !> @brief Unstructured grid discretization
  type, extends(disbasetype) :: disutype
    integer(I4B), pointer :: njausr => null() ! user-specified nja size
    integer(I4B), pointer :: nvert => null() ! number of x,y vertices
    integer(I4B), pointer :: nangldegxerr => null() ! number of horizontal connections with an ANGLDEGX value that is inconsistent with the reverse connection or with the vertices
    real(dp), pointer :: voffsettol => null() ! vertical offset tolerance
    real(dp), dimension(:, :), pointer, contiguous :: vertices => null() ! cell vertices stored as 2d array of x and y
    real(dp), dimension(:, :), pointer, contiguous :: cellxy => null() ! cell center stored as 2d array of x and y
    real(dp), dimension(:), pointer, contiguous :: top1d => null() ! (size:nodesuser) cell top elevation
    real(dp), dimension(:), pointer, contiguous :: bot1d => null() ! (size:nodesuser) cell bottom elevation
    real(dp), dimension(:), pointer, contiguous :: area1d => null() ! (size:nodesuser) cell area, in plan view
    integer(I4B), dimension(:), pointer, contiguous :: iainp => null() ! (size:nodesuser+1) user iac converted ia
    integer(I4B), dimension(:), pointer, contiguous :: jainp => null() ! (size:njausr) user-input ja array
    integer(I4B), dimension(:), pointer, contiguous :: ihcinp => null() ! (size:njausr) user-input ihc array
    real(dp), dimension(:), pointer, contiguous :: cl12inp => null() ! (size:njausr) user-input cl12 array
    real(dp), dimension(:), pointer, contiguous :: hwvainp => null() ! (size:njausr) user-input hwva array
    real(dp), dimension(:), pointer, contiguous :: angldegxinp => null() ! (size:njausr) user-input angldegx array
    integer(I4B), pointer :: iangledegx => null() ! =1 when angle information was present in input, 0 otherwise
    integer(I4B), dimension(:), pointer, contiguous :: iavert => null() ! cell vertex pointer ia array
    integer(I4B), dimension(:), pointer, contiguous :: javert => null() ! cell vertex pointer ja array
    integer(I4B), dimension(:), pointer, contiguous :: idomain => null() ! idomain (nodes)
    logical(LGP) :: readfromfile ! True, when DIS is read from file (almost always)
 
  contains
 
    procedure :: dis_df => disu_df
    procedure :: disu_load
    procedure :: dis_da => disu_da
    procedure :: get_dis_type => get_dis_type
    procedure :: get_dis_enum => get_dis_enum
    procedure :: disu_ck
    procedure :: disu_write_warning
    procedure :: disu_face_normal
    procedure :: grid_finalize
    procedure :: get_nodenumber_idx1
    procedure :: nodeu_to_string
    procedure :: nodeu_to_array
    procedure :: nodeu_from_string
    procedure :: nodeu_from_cellid
    procedure :: connection_normal
    procedure :: connection_vector
    procedure :: supports_layers
    procedure :: get_ncpl
    procedure :: get_polyverts
    procedure :: get_npolyverts
    procedure :: get_max_npolyverts
    procedure, public :: record_array
    procedure, public :: record_srcdst_list_header
    ! -- private
    procedure :: allocate_scalars
    procedure :: allocate_arrays
    procedure :: allocate_arrays_mem
    procedure :: source_options
    procedure :: source_dimensions
    procedure :: source_griddata
    procedure :: source_connectivity
    procedure :: source_vertices
    procedure :: source_cell2d
    procedure :: log_options
    procedure :: log_dimensions
    procedure :: log_griddata
    procedure :: log_connectivity
    procedure :: define_cellverts
    procedure :: write_grb
    !
    ! -- Read a node-sized model array (reduced or not)
    procedure :: read_int_array
    procedure :: read_dbl_array
 
  end type disutype
 
  type disufoundtype
    logical :: length_units = .false.
    logical :: nogrb = .false.
    logical :: xorigin = .false.
    logical :: yorigin = .false.
    logical :: angrot = .false.
    logical :: voffsettol = .false.
    logical :: nodes = .false.
    logical :: nja = .false.
    logical :: nvert = .false.
    logical :: top = .false.
    logical :: bot = .false.
    logical :: area = .false.
    logical :: idomain = .false.
    logical :: iac = .false.
    logical :: ja = .false.
    logical :: ihc = .false.
    logical :: cl12 = .false.
    logical :: hwva = .false.
    logical :: angldegx = .false.
    logical :: iv = .false.
    logical :: xv = .false.
    logical :: yv = .false.
    logical :: icell2d = .false.
    logical :: xc = .false.
    logical :: yc = .false.
    logical :: ncvert = .false.
    logical :: icvert = .false.
  end type disufoundtype
 
contains
 
  !> @brief Create a new unstructured discretization object
  !<
  subroutine disu_cr(dis, name_model, input_mempath, inunit, iout)
    ! -- dummy
    class(disbasetype), pointer :: dis
    character(len=*), intent(in) :: name_model
    character(len=*), intent(in) :: input_mempath
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    ! -- local
    type(disutype), pointer :: disnew
    character(len=*), parameter :: fmtheader = &
      "(1X, /1X, 'DISU -- UNSTRUCTURED GRID DISCRETIZATION PACKAGE,', &
      &' VERSION 2 : 3/27/2014 - INPUT READ FROM MEMPATH: ', A, //)"
    !
    ! -- Create a new discretization object
    allocate (disnew)
    dis => disnew
    !
    ! -- Allocate scalars and assign data
    call dis%allocate_scalars(name_model, input_mempath)
    dis%inunit = inunit
    dis%iout = iout
    !
    ! -- If disu is enabled
    if (inunit > 0) then
      !
      ! -- Identify package
      if (iout > 0) then
        write (iout, fmtheader) dis%input_mempath
      end if
      !
      ! -- load disu
      call disnew%disu_load()
    end if
    !
  end subroutine disu_cr
 
  !> @brief Transfer IDM data into this discretization object
  !<
  subroutine disu_load(this)
    ! -- dummy
    class(disutype) :: this
    !
    ! -- source input data
    call this%source_options()
    call this%source_dimensions()
    call this%source_griddata()
    call this%source_connectivity()
    !
    ! -- If NVERT specified and greater than 0, then source VERTICES and CELL2D
    if (this%nvert > 0) then
      call this%source_vertices()
      call this%source_cell2d()
    else
      ! -- connection direction information cannot be calculated
      this%icondir = 0
    end if
    !
    ! -- Make some final disu checks on the non-reduced user-provided
    !    input
    call this%disu_ck()
    !
  end subroutine disu_load
 
  !> @brief Define the discretization
  !<
  subroutine disu_df(this)
    ! -- dummy
    class(disutype) :: this
    !
    call this%grid_finalize()
    !
  end subroutine disu_df
 
  !> @brief Finalize the grid
  !<
  subroutine grid_finalize(this)
    ! -- dummy
    class(disutype) :: this
    ! -- locals
    integer(I4B) :: n
    integer(I4B) :: node
    integer(I4B) :: noder
    integer(I4B) :: nrsize
    ! -- formats
    character(len=*), parameter :: fmtdz = &
      "('CELL (',i0,',',i0,',',i0,') THICKNESS <= 0. ', &
      &'TOP, BOT: ',2(1pg24.15))"
    character(len=*), parameter :: fmtnr = &
      "(/1x, 'The specified IDOMAIN results in a reduced number of cells.',&
      &/1x, 'Number of user nodes: ',I0,&
      &/1X, 'Number of nodes in solution: ', I0, //)"
    !
    ! -- count active cells
    this%nodes = 0
    do n = 1, this%nodesuser
      if (this%idomain(n) > 0) this%nodes = this%nodes + 1
    end do
    !
    ! -- Check to make sure nodes is a valid number
    if (this%nodes == 0) then
      call store_error('Model does not have any active nodes. &
                       &Ensure IDOMAIN array has some values greater &
                       &than zero.')
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- Write message if reduced grid
    if (this%nodes < this%nodesuser) then
      write (this%iout, fmtnr) this%nodesuser, this%nodes
    end if
    !
    ! -- Array size is now known, so allocate
    call this%allocate_arrays()
    !
    ! -- Fill the nodereduced array with the reduced nodenumber, or
    !    a negative number to indicate it is a pass-through cell, or
    !    a zero to indicate that the cell is excluded from the
    !    solution. (negative idomain not supported for disu)
    if (this%nodes < this%nodesuser) then
      noder = 1
      do node = 1, this%nodesuser
        if (this%idomain(node) > 0) then
          this%nodereduced(node) = noder
          noder = noder + 1
        elseif (this%idomain(node) < 0) then
          this%nodereduced(node) = -1
        else
          this%nodereduced(node) = 0
        end if
      end do
    end if
    !
    ! -- Fill nodeuser if a reduced grid
    if (this%nodes < this%nodesuser) then
      noder = 1
      do node = 1, this%nodesuser
        if (this%idomain(node) > 0) then
          this%nodeuser(noder) = node
          noder = noder + 1
        end if
      end do
    end if
    !
    ! -- Move top1d, bot1d, and area1d into top, bot, and area
    do node = 1, this%nodesuser
      noder = node
      if (this%nodes < this%nodesuser) noder = this%nodereduced(node)
      if (noder <= 0) cycle
      this%top(noder) = this%top1d(node)
      this%bot(noder) = this%bot1d(node)
      this%area(noder) = this%area1d(node)
    end do
    !
    ! -- fill cell center coordinates
    if (this%nvert > 0) then
      do node = 1, this%nodesuser
        noder = node
        if (this%nodes < this%nodesuser) noder = this%nodereduced(node)
        if (noder <= 0) cycle
        this%xc(noder) = this%cellxy(1, node)
        this%yc(noder) = this%cellxy(2, node)
      end do
    else
      call mem_reallocate(this%xc, 0, 'XC', this%memoryPath)
      call mem_reallocate(this%yc, 0, 'YC', this%memoryPath)
    end if
    !
    ! -- create and fill the connections object
    nrsize = 0
    if (this%nodes < this%nodesuser) nrsize = this%nodes
    allocate (this%con)
    call this%con%disuconnections(this%name_model, this%nodes, &
                                  this%nodesuser, nrsize, &
                                  this%nodereduced, this%nodeuser, &
                                  this%iainp, this%jainp, &
                                  this%ihcinp, this%cl12inp, &
                                  this%hwvainp, this%angldegxinp, &
                                  this%iangledegx)
    this%nja = this%con%nja
    this%njas = this%con%njas
    !
  end subroutine grid_finalize
 
  !> @brief Check discretization info
  !<
  subroutine disu_ck(this)
    ! -- dummy
    class(disutype) :: this
    ! -- local
    integer(I4B) :: n, m
    integer(I4B) :: ipos, jpos, kpos
    integer(I4B) :: ihc
    integer(I4B) :: nsym, ndir, nvrt, nwarn
    logical :: found
    real(DP) :: dz
    real(DP) :: angn, angm, dang
    real(DP) :: dx, dy, dist, cosang, angc
    real(DP) :: xnorm, ynorm, angv
    real(DP) :: angn_sym, angm_sym
    real(DP) :: angn_dir, angc_dir
    real(DP) :: angn_vrt, angv_vrt
    real(DP) :: angn_warn, angc_warn
    integer(I4B) :: n_sym, m_sym
    integer(I4B) :: n_dir, m_dir
    integer(I4B) :: n_vrt, m_vrt
    integer(I4B) :: n_warn, m_warn
    character(len=LENBIGLINE) :: bigmsg
    real(DP), parameter :: angtol = 0.01_dp !< tolerance (degrees) for the 180-degree reciprocity check
    real(DP), parameter :: angvtol = 1.0_dp !< tolerance (degrees) for the comparison with the normal computed from the vertices
    real(DP), parameter :: cos45 = 0.70710678118654752_dp
    ! -- formats
    character(len=*), parameter :: fmtidm = &
      &"('Invalid idomain value ', i0, ' specified for node ', i0)"
    character(len=*), parameter :: fmtangsym = &
      "('ANGLDEGX values for ', i0, ' cell faces associated with &
      &horizontal connections in the DISU Package are inconsistent with &
      &the values for the reverse connections (for example, ANGLDEGX = ', &
      &f0.3, ' for the connection from cell ', i0, ' to cell ', i0, &
      &' and ANGLDEGX = ', f0.3, ' for the connection from cell ', i0, &
      &' to cell ', i0, '). The two values must differ by 180 degrees &
      &(within 0.01 degrees) because they are the outward normals of the &
      &two sides of the same cell face. Only the value for the connection &
      &from the lower to the higher cell number is used.')"
    character(len=*), parameter :: fmtangdir = &
      "('ANGLDEGX values for ', i0, ' cell faces associated with &
      &horizontal connections in the DISU Package point away from the &
      &connected cell (for example, ANGLDEGX = ', f0.3, ' for the &
      &connection from cell ', i0, ' to cell ', i0, ', but the direction &
      &from the center of cell ', i0, ' to the center of cell ', i0, &
      &' is ', f0.3, ' degrees). The centers of two connected cells must &
      &lie on opposite sides of their shared face, so ANGLDEGX, the &
      &outward normal of the face, must be within 90 degrees of the &
      &direction from the cell center to the center of the connected &
      &cell, in the coordinate system of the VERTICES.')"
    character(len=*), parameter :: fmtangvrt = &
      "('ANGLDEGX values for ', i0, ' cell faces associated with &
      &horizontal connections in the DISU Package differ by more than 1 &
      &degree from the outward normal of the face computed from the two &
      &VERTICES shared by the connected cells (for example, ANGLDEGX = ', &
      &f0.3, ' for the connection from cell ', i0, ' to cell ', i0, &
      &', but the normal computed from the vertices is ', f0.3, &
      &' degrees). ANGLDEGX must be expressed in the coordinate system of &
      &the VERTICES.')"
    character(len=*), parameter :: fmtangwarn = &
      "('ANGLDEGX values for ', i0, ' cell faces associated with &
      &horizontal connections in the DISU Package deviate by more than 45 &
      &degrees from the direction between the centers of the connected &
      &cells (for example, ANGLDEGX = ', f0.3, ' for the connection from &
      &cell ', i0, ' to cell ', i0, ', but the direction between the cell &
      &centers is ', f0.3, ' degrees). This is not necessarily an error, &
      &but such connections depart substantially from the assumptions of &
      &the control-volume finite-difference method and may warrant a &
      &closer look at the grid.')"
    character(len=*), parameter :: fmtdz = &
      &"('Cell ', i0, ' with thickness <= 0. Top, bot: ', 2(1pg24.15))"
    character(len=*), parameter :: fmtarea = &
      &"('Cell ', i0, ' with area <= 0. Area: ', 1(1pg24.15))"
    character(len=*), parameter :: fmtjan = &
      &"('Cell ', i0, ' must have its first connection be itself.  Found: ', i0)"
    character(len=*), parameter :: fmtjam = &
      &"('Cell ', i0, ' has invalid connection in JA.  Found: ', i0)"
    character(len=*), parameter :: fmterrmsg = &
      "('Top elevation (', 1pg15.6, ') for cell ', i0, ' is above bottom &
      &elevation (', 1pg15.6, ') for cell ', i0, '. Based on node numbering &
      &rules cell ', i0, ' must be below cell ', i0, '.')"
    !
    ! -- Check connectivity
    do n = 1, this%nodesuser
      !
      ! -- Ensure first connection is to itself, and
      !    that ja(ia(n)) is positive
      ipos = this%iainp(n)
      m = this%jainp(ipos)
      if (m < 0) then
        m = abs(m)
        this%jainp(ipos) = m
      end if
      if (n /= m) then
        write (errmsg, fmtjan) n, m
        call store_error(errmsg)
      end if
      !
      ! -- Check for valid node numbers in connected cells
      do ipos = this%iainp(n) + 1, this%iainp(n + 1) - 1
        m = this%jainp(ipos)
        if (m < 0 .or. m > this%nodesuser) then
          ! -- make sure first connection is to itself
          write (errmsg, fmtjam) n, m
          call store_error(errmsg)
        end if
      end do
    end do
    !
    ! -- terminate if errors found
    if (count_errors() > 0) then
      if (this%inunit > 0) then
        call store_error_filename(this%input_fname)
      end if
    end if
    !
    ! -- Ensure idomain values are valid
    do n = 1, this%nodesuser
      if (this%idomain(n) > 1 .or. this%idomain(n) < 0) then
        write (errmsg, fmtidm) this%idomain(n), n
        call store_error(errmsg)
      end if
    end do
    !
    ! -- Check for zero and negative thickness and zero or negative areas
    !    for cells with idomain == 1
    do n = 1, this%nodesuser
      if (this%idomain(n) == 1) then
        dz = this%top1d(n) - this%bot1d(n)
        if (dz <= dzero) then
          write (errmsg, fmt=fmtdz) n, this%top1d(n), this%bot1d(n)
          call store_error(errmsg)
        end if
        if (this%area1d(n) <= dzero) then
          write (errmsg, fmt=fmtarea) n, this%area1d(n)
          call store_error(errmsg)
        end if
      end if
    end do
    !
    ! -- check to make sure voffsettol is >= 0
    if (this%voffsettol < dzero) then
      write (errmsg, '(a, 1pg15.6)') &
        'Vertical offset tolerance must be greater than zero. Found ', &
        this%voffsettol
      call store_error(errmsg)
      if (this%inunit > 0) then
        call store_error_filename(this%input_fname)
      end if
    end if
    !
    ! -- For cell n, ensure that underlying cells have tops less than
    !    or equal to the bottom of cell n
    do n = 1, this%nodesuser
      do ipos = this%iainp(n) + 1, this%iainp(n + 1) - 1
        m = this%jainp(ipos)
        ihc = this%ihcinp(ipos)
        if (ihc == 0 .and. m > n) then
          dz = this%top1d(m) - this%bot1d(n)
          if (dz > this%voffsettol) then
            write (errmsg, fmterrmsg) this%top1d(m), m, this%bot1d(n), n, m, n
            call store_error(errmsg)
          end if
        end if
      end do
    end do
    !
    ! -- Check ANGLDEGX for horizontal connections between active cells.
    !    ANGLDEGX is the outward normal of the shared face, so the values
    !    for a connection and its reverse connection must differ by 180
    !    degrees.  If VERTICES and CELL2D are available, the normal must
    !    also point into the connected cell (the centers of the two cells
    !    lie on opposite sides of the shared face) and must agree with the
    !    normal computed from the two vertices shared by the cells.  These
    !    three conditions are reported as warnings here and counted in
    !    nangldegxerr; the NPF Package terminates with an error if the count
    !    is nonzero and ANGLDEGX is required input (XT3D, K22, or
    !    SAVE_SPECIFIC_DISCHARGE), because ANGLDEGX has no effect otherwise.
    !    A fourth check flags normals that deviate by more than 45 degrees
    !    from the direction between the cell centers.  That is not an error,
    !    but such connections depart from the assumptions of the CVFD
    !    method and are reported as a warning only.
    if (this%iangledegx == 1) then
      nsym = 0
      ndir = 0
      nvrt = 0
      nwarn = 0
      n_sym = 0
      m_sym = 0
      n_dir = 0
      m_dir = 0
      n_vrt = 0
      m_vrt = 0
      n_warn = 0
      m_warn = 0
      angn_sym = dzero
      angm_sym = dzero
      angn_dir = dzero
      angc_dir = dzero
      angn_vrt = dzero
      angv_vrt = dzero
      angn_warn = dzero
      angc_warn = dzero
      do n = 1, this%nodesuser
        if (this%idomain(n) == 0) cycle
        do ipos = this%iainp(n) + 1, this%iainp(n + 1) - 1
          m = this%jainp(ipos)
          if (m < 1 .or. m > this%nodesuser) cycle
          if (this%ihcinp(ipos) == 0) cycle
          if (this%idomain(m) == 0) cycle
          angn = this%angldegxinp(ipos)
          !
          ! -- the reverse connection must have the opposite normal;
          !    check each pair once
          if (m > n) then
            jpos = 0
            do kpos = this%iainp(m) + 1, this%iainp(m + 1) - 1
              if (this%jainp(kpos) == n) then
                jpos = kpos
                exit
              end if
            end do
            if (jpos > 0) then
              angm = this%angldegxinp(jpos)
              dang = modulo(angm - angn, 360.0_dp)
              if (abs(dang - 180.0_dp) > angtol) then
                nsym = nsym + 1
                if (nsym == 1) then
                  n_sym = n
                  m_sym = m
                  angn_sym = angn
                  angm_sym = angm
                end if
              end if
            end if
          end if
          !
          ! -- the normal must point from cell n toward cell m
          if (this%nvert > 0) then
            dx = this%cellxy(1, m) - this%cellxy(1, n)
            dy = this%cellxy(2, m) - this%cellxy(2, n)
            dist = sqrt(dx * dx + dy * dy)
            if (dist > dzero) then
              cosang = (cos(angn * dpio180) * dx + &
                        sin(angn * dpio180) * dy) / dist
              angc = modulo(atan2(dy, dx) / dpio180, 360.0_dp)
              if (cosang <= dzero) then
                ndir = ndir + 1
                if (ndir == 1) then
                  n_dir = n
                  m_dir = m
                  angn_dir = angn
                  angc_dir = angc
                end if
              else if (cosang < cos45) then
                nwarn = nwarn + 1
                if (nwarn == 1) then
                  angn_warn = angn
                  angc_warn = angc
                  n_warn = n
                  m_warn = m
                end if
              end if
            end if
            !
            ! -- the normal must agree with the normal of the face computed
            !    from the vertices; check each pair once, and skip pairs
            !    that do not share exactly two vertices
            if (m > n) then
              call this%disu_face_normal(n, m, found, xnorm, ynorm)
              if (found) then
                angv = modulo(atan2(ynorm, xnorm) / dpio180, 360.0_dp)
                dang = abs(modulo(angn - angv + 180.0_dp, 360.0_dp) - 180.0_dp)
                if (dang > angvtol) then
                  nvrt = nvrt + 1
                  if (nvrt == 1) then
                    n_vrt = n
                    m_vrt = m
                    angn_vrt = angn
                    angv_vrt = angv
                  end if
                end if
              end if
            end if
          end if
        end do
      end do
      this%nangldegxerr = nsym + ndir + nvrt
      !
      ! -- store warnings for the final report and also write them to the
      !    model listing file now, so that they are available if the run
      !    later terminates abnormally because of an invalid matrix
      if (nsym > 0) then
        write (bigmsg, fmtangsym) nsym, angn_sym, n_sym, m_sym, &
          angm_sym, m_sym, n_sym
        call store_warning(bigmsg)
        call this%disu_write_warning(bigmsg)
      end if
      if (ndir > 0) then
        write (bigmsg, fmtangdir) ndir, angn_dir, n_dir, m_dir, &
          n_dir, m_dir, angc_dir
        call store_warning(bigmsg)
        call this%disu_write_warning(bigmsg)
      end if
      if (nvrt > 0) then
        write (bigmsg, fmtangvrt) nvrt, angn_vrt, n_vrt, m_vrt, angv_vrt
        call store_warning(bigmsg)
        call this%disu_write_warning(bigmsg)
      end if
      if (nwarn > 0) then
        write (bigmsg, fmtangwarn) nwarn, angn_warn, n_warn, m_warn, angc_warn
        call store_warning(bigmsg)
        call this%disu_write_warning(bigmsg)
      end if
    end if
    !
    ! -- terminate if errors found
    if (count_errors() > 0) then
      if (this%inunit > 0) then
        call store_error_filename(this%input_fname)
      end if
    end if
    !
  end subroutine disu_ck
 
  !> @brief Write a warning message to the model listing file
  !<
  subroutine disu_write_warning(this, msg)
    ! -- dummy
    class(disutype) :: this
    character(len=*), intent(in) :: msg
    !
    if (this%iout > 0) then
      call write_message_counter('WARNING: '//trim(msg), iunit=this%iout, &
                                 skipbefore=1, skipafter=1)
    end if
  end subroutine disu_write_warning
 
  !> @brief Outward normal of the face shared by two cells from the vertices
  !!
  !! Finds the vertices shared by cells n and m.  Vertices are matched by
  !! coordinates rather than by vertex number, so that grids that list a
  !! separate set of vertices for every cell are handled the same as grids
  !! with shared vertex numbers.  If exactly two vertices are shared, they
  !! define the shared face, and the unit normal to that face that points
  !! from cell n toward cell m is returned with found set to true.
  !! Otherwise found is set to false.  Requires VERTICES and CELL2D.
  !<
  subroutine disu_face_normal(this, n, m, found, xnorm, ynorm)
    ! -- dummy
    class(disutype) :: this
    integer(I4B), intent(in) :: n !< cell (user node number)
    integer(I4B), intent(in) :: m !< connected cell (user node number)
    logical, intent(out) :: found !< true if the cells share exactly two vertices
    real(DP), intent(out) :: xnorm !< x component of the unit normal from n toward m
    real(DP), intent(out) :: ynorm !< y component of the unit normal from n toward m
    ! -- local
    integer(I4B) :: i, j, iv, jv, nshared
    integer(I4B), dimension(2) :: ivshared
    real(DP) :: ex, ey, elen, dx, dy, tol
    !
    found = .false.
    xnorm = dzero
    ynorm = dzero
    nshared = 0
    !
    ! -- Vertices are shared if their coordinates agree to within a small
    !    fraction of the distance between the two cell centers
    dx = this%cellxy(1, m) - this%cellxy(1, n)
    dy = this%cellxy(2, m) - this%cellxy(2, n)
    tol = dem6 * sqrt(dx * dx + dy * dy)
    !
    ! -- Count the distinct vertices of cell n that also belong to cell m
    !    (the vertex list of a cell repeats its first vertex at the end)
    do i = this%iavert(n), this%iavert(n + 1) - 1
      iv = this%javert(i)
      if (any(this%javert(this%iavert(n):i - 1) == iv)) cycle
      do j = this%iavert(m), this%iavert(m + 1) - 1
        jv = this%javert(j)
        if (abs(this%vertices(1, jv) - this%vertices(1, iv)) <= tol .and. &
            abs(this%vertices(2, jv) - this%vertices(2, iv)) <= tol) then
          nshared = nshared + 1
          if (nshared > 2) return
          ivshared(nshared) = iv
          exit
        end if
      end do
    end do
    if (nshared /= 2) return
    !
    ! -- Unit normal to the shared face, oriented from cell n toward cell m
    ex = this%vertices(1, ivshared(2)) - this%vertices(1, ivshared(1))
    ey = this%vertices(2, ivshared(2)) - this%vertices(2, ivshared(1))
    elen = sqrt(ex * ex + ey * ey)
    if (elen <= dzero) return
    xnorm = ey / elen
    ynorm = -ex / elen
    if (xnorm * dx + ynorm * dy < dzero) then
      xnorm = -xnorm
      ynorm = -ynorm
    end if
    found = .true.
  end subroutine disu_face_normal
 
  !> @brief Deallocate variables
  !<
  subroutine disu_da(this)
    ! -- dummy
    class(disutype) :: this
    !
    ! -- Deallocate idm memory
    call memorystore_remove(this%name_model, 'DISU', idm_context)
    call memorystore_remove(component=this%name_model, &
                            context=idm_context)
    !
    ! -- scalars
    call mem_deallocate(this%njausr)
    call mem_deallocate(this%nvert)
    call mem_deallocate(this%nangldegxerr)
    call mem_deallocate(this%voffsettol)
    call mem_deallocate(this%iangledegx)
    !
    ! -- arrays
    if (this%readFromFile) then
      call mem_deallocate(this%top1d)
      call mem_deallocate(this%bot1d)
      call mem_deallocate(this%area1d)
      if (associated(this%iavert)) then
        call mem_deallocate(this%iavert)
        call mem_deallocate(this%javert)
      end if
      call mem_deallocate(this%vertices)
      call mem_deallocate(this%iainp)
      call mem_deallocate(this%jainp)
      call mem_deallocate(this%ihcinp)
      call mem_deallocate(this%cl12inp)
      call mem_deallocate(this%hwvainp)
      call mem_deallocate(this%angldegxinp)
    end if
    !
    call mem_deallocate(this%idomain)
    call mem_deallocate(this%cellxy)
    !
    call mem_deallocate(this%nodeuser)
    call mem_deallocate(this%nodereduced)
    !
    ! -- DisBaseType deallocate
    call this%DisBaseType%dis_da()
    !
  end subroutine disu_da
 
  !> @brief Convert a user nodenumber to a string (nodenumber)
  !<
  subroutine nodeu_to_string(this, nodeu, str)
    ! -- dummy
    class(disutype) :: this
    integer(I4B), intent(in) :: nodeu
    character(len=*), intent(inout) :: str
    ! -- local
    character(len=10) :: nstr
    !
    write (nstr, '(i0)') nodeu
    str = '('//trim(adjustl(nstr))//')'
    !
  end subroutine nodeu_to_string
 
  !> @brief Convert a user nodenumber to an array (nodenumber)
  !<
  subroutine nodeu_to_array(this, nodeu, arr)
    class(disutype) :: this
    integer(I4B), intent(in) :: nodeu
    integer(I4B), dimension(:), intent(inout) :: arr
    ! -- local
    integer(I4B) :: isize
    !
    ! -- check the size of arr
    isize = size(arr)
    if (isize /= this%ndim) then
      write (errmsg, '(a,i0,a,i0,a)') &
        'Program error: nodeu_to_array size of array (', isize, &
        ') is not equal to the discretization dimension (', this%ndim, ')'
      call store_error(errmsg, terminate=.true.)
    end if
    !
    ! -- fill array
    arr(1) = nodeu
    !
  end subroutine nodeu_to_array
 
  !> @brief Copy options from IDM into package
  !<
  subroutine source_options(this)
    ! -- dummy
    class(disutype) :: this
    ! -- locals
    character(len=LENVARNAME), dimension(3) :: lenunits = &
      &[character(len=LENVARNAME) :: 'FEET', 'METERS', 'CENTIMETERS']
    type(disufoundtype) :: found
    !
    ! -- update defaults with idm sourced values
    call mem_set_value(this%lenuni, 'LENGTH_UNITS', this%input_mempath, &
                       lenunits, found%length_units)
    call mem_set_value(this%nogrb, 'NOGRB', this%input_mempath, found%nogrb)
    call mem_set_value(this%xorigin, 'XORIGIN', this%input_mempath, found%xorigin)
    call mem_set_value(this%yorigin, 'YORIGIN', this%input_mempath, found%yorigin)
    call mem_set_value(this%angrot, 'ANGROT', this%input_mempath, found%angrot)
    call mem_set_value(this%voffsettol, 'VOFFSETTOL', this%input_mempath, &
                       found%voffsettol)
    !
    ! -- log values to list file
    if (this%iout > 0) then
      call this%log_options(found)
    end if
    !
  end subroutine source_options
 
  !> @brief Write user options to list file
  !<
  subroutine log_options(this, found)
    ! -- dummy
    class(disutype) :: this
    type(disufoundtype), intent(in) :: found
    !
    write (this%iout, '(1x,a)') 'Setting Discretization Options'
    !
    if (found%length_units) then
      write (this%iout, '(4x,a,i0)') 'Model length unit [0=UND, 1=FEET, &
      &2=METERS, 3=CENTIMETERS] set as ', this%lenuni
    end if
    !
    if (found%nogrb) then
      write (this%iout, '(4x,a,i0)') 'Binary grid file [0=GRB, 1=NOGRB] &
        &set as ', this%nogrb
    end if
    !
    if (found%xorigin) then
      write (this%iout, '(4x,a,G0)') 'XORIGIN = ', this%xorigin
    end if
    !
    if (found%yorigin) then
      write (this%iout, '(4x,a,G0)') 'YORIGIN = ', this%yorigin
    end if
    !
    if (found%angrot) then
      write (this%iout, '(4x,a,G0)') 'ANGROT = ', this%angrot
    end if
    !
    if (found%voffsettol) then
      write (this%iout, '(4x,a,G0)') 'VERTICAL_OFFSET_TOLERANCE = ', &
        this%voffsettol
    end if
    !
    write (this%iout, '(1x,a,/)') 'End Setting Discretization Options'
    !
  end subroutine log_options
 
  !> @brief Copy dimensions from IDM into package
  !<
  subroutine source_dimensions(this)
    ! -- dummy
    class(disutype) :: this
    ! -- locals
    integer(I4B) :: n
    type(disufoundtype) :: found
    !
    ! -- update defaults with idm sourced values
    call mem_set_value(this%nodesuser, 'NODES', this%input_mempath, found%nodes)
    call mem_set_value(this%njausr, 'NJA', this%input_mempath, found%nja)
    call mem_set_value(this%nvert, 'NVERT', this%input_mempath, found%nvert)
    !
    ! -- log simulation values
    if (this%iout > 0) then
      call this%log_dimensions(found)
    end if
    !
    ! -- verify dimensions were set
    if (this%nodesuser < 1) then
      call store_error( &
        'NODES was not specified or was specified incorrectly.')
    end if
    if (this%njausr < 1) then
      call store_error( &
        'NJA was not specified or was specified incorrectly.')
    end if
    !
    ! -- terminate if errors were detected
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- allocate vectors that are the size of nodesuser
    this%readFromFile = .true.
    call mem_allocate(this%top1d, this%nodesuser, 'TOP1D', this%memoryPath)
    call mem_allocate(this%bot1d, this%nodesuser, 'BOT1D', this%memoryPath)
    call mem_allocate(this%area1d, this%nodesuser, 'AREA1D', this%memoryPath)
    call mem_allocate(this%idomain, this%nodesuser, 'IDOMAIN', this%memoryPath)
    call mem_allocate(this%vertices, 2, this%nvert, 'VERTICES', this%memoryPath)
    call mem_allocate(this%iainp, this%nodesuser + 1, 'IAINP', this%memoryPath)
    call mem_allocate(this%jainp, this%njausr, 'JAINP', this%memoryPath)
    call mem_allocate(this%ihcinp, this%njausr, 'IHCINP', this%memoryPath)
    call mem_allocate(this%cl12inp, this%njausr, 'CL12INP', this%memoryPath)
    call mem_allocate(this%hwvainp, this%njausr, 'HWVAINP', this%memoryPath)
    call mem_allocate(this%angldegxinp, this%njausr, 'ANGLDEGXINP', &
                      this%memoryPath)
    if (this%nvert > 0) then
      call mem_allocate(this%cellxy, 2, this%nodesuser, 'CELLXY', this%memoryPath)
    else
      call mem_allocate(this%cellxy, 2, 0, 'CELLXY', this%memoryPath)
    end if
    !
    ! -- initialize all cells to be active (idomain = 1)
    do n = 1, this%nodesuser
      this%idomain(n) = 1
    end do
    !
  end subroutine source_dimensions
 
  !> @brief Write dimensions to list file
  !<
  subroutine log_dimensions(this, found)
    class(disutype) :: this
    type(disufoundtype), intent(in) :: found
    !
    write (this%iout, '(1x,a)') 'Setting Discretization Dimensions'
    !
    if (found%nodes) then
      write (this%iout, '(4x,a,i0)') 'NODES = ', this%nodesuser
    end if
    !
    if (found%nja) then
      write (this%iout, '(4x,a,i0)') 'NJA = ', this%njausr
    end if
    !
    if (found%nvert) then
      write (this%iout, '(4x,a,i0)') 'NVERT = ', this%nvert
    end if
    !
    write (this%iout, '(1x,a,/)') 'End Setting Discretization Dimensions'
    !
  end subroutine log_dimensions
 
  !> @brief Copy grid data from IDM into package
  !<
  subroutine source_griddata(this)
    ! -- dummy
    class(disutype) :: this
    ! -- locals
    type(disufoundtype) :: found
    !
    ! -- update defaults with idm sourced values
    call mem_set_value(this%top1d, 'TOP', this%input_mempath, found%top)
    call mem_set_value(this%bot1d, 'BOT', this%input_mempath, found%bot)
    call mem_set_value(this%area1d, 'AREA', this%input_mempath, found%area)
    call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain)
    !
    ! -- log simulation values
    if (this%iout > 0) then
      call this%log_griddata(found)
    end if
    !
  end subroutine source_griddata
 
  !> @brief Write griddata found to list file
  !<
  subroutine log_griddata(this, found)
    ! -- dummy
    class(disutype) :: this
    type(disufoundtype), intent(in) :: found
    !
    write (this%iout, '(1x,a)') 'Setting Discretization Griddata'
    !
    if (found%top) then
      write (this%iout, '(4x,a)') 'TOP set from input file'
    end if
    !
    if (found%bot) then
      write (this%iout, '(4x,a)') 'BOT set from input file'
    end if
    !
    if (found%area) then
      write (this%iout, '(4x,a)') 'AREA set from input file'
    end if
    !
    if (found%idomain) then
      write (this%iout, '(4x,a)') 'IDOMAIN set from input file'
    end if
    !
    write (this%iout, '(1x,a,/)') 'End Setting Discretization Griddata'
    !
  end subroutine log_griddata
 
  !> @brief Copy grid connectivity info from IDM into package
  !<
  subroutine source_connectivity(this)
    ! -- dummy
    class(disutype) :: this
    ! -- locals
    type(disufoundtype) :: found
    integer(I4B), dimension(:), contiguous, pointer :: iac => null()
    ! -- formats
    !
    ! -- update defaults with idm sourced values
    call mem_set_value(this%jainp, 'JA', this%input_mempath, found%ja)
    call mem_set_value(this%ihcinp, 'IHC', this%input_mempath, found%ihc)
    call mem_set_value(this%cl12inp, 'CL12', this%input_mempath, found%cl12)
    call mem_set_value(this%hwvainp, 'HWVA', this%input_mempath, found%hwva)
    call mem_set_value(this%angldegxinp, 'ANGLDEGX', this%input_mempath, &
                       found%angldegx)
    !
    ! -- set pointer to iac input array
    call mem_setptr(iac, 'IAC', this%input_mempath)
    !
    ! -- Convert iac to ia
    if (associated(iac)) call iac_to_ia(iac, this%iainp)
    !
    ! -- Set angldegx flag if found
    if (found%angldegx) this%iangledegx = 1
    !
    ! -- log simulation values
    if (this%iout > 0) then
      call this%log_connectivity(found, iac)
    end if
    !
  end subroutine source_connectivity
 
  !> @brief Write griddata found to list file
  !<
  subroutine log_connectivity(this, found, iac)
    class(disutype) :: this
    type(disufoundtype), intent(in) :: found
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: iac
    !
    write (this%iout, '(1x,a)') 'Setting Discretization Connectivity'
    !
    if (associated(iac)) then
      write (this%iout, '(4x,a)') 'IAC set from input file'
    end if
    !
    if (found%ja) then
      write (this%iout, '(4x,a)') 'JA set from input file'
    end if
    !
    if (found%ihc) then
      write (this%iout, '(4x,a)') 'IHC set from input file'
    end if
    !
    if (found%cl12) then
      write (this%iout, '(4x,a)') 'CL12 set from input file'
    end if
    !
    if (found%hwva) then
      write (this%iout, '(4x,a)') 'HWVA set from input file'
    end if
    !
    if (found%angldegx) then
      write (this%iout, '(4x,a)') 'ANGLDEGX set from input file'
    end if
    !
    write (this%iout, '(1x,a,/)') 'End Setting Discretization Connectivity'
    !
  end subroutine log_connectivity
 
  !> @brief Copy grid vertex data from IDM into package
  !<
  subroutine source_vertices(this)
    ! -- dummy
    class(disutype) :: this
    ! -- local
    integer(I4B) :: i
    real(DP), dimension(:), contiguous, pointer :: vert_x => null()
    real(DP), dimension(:), contiguous, pointer :: vert_y => null()
    ! -- formats
    !
    ! -- set pointers to memory manager input arrays
    call mem_setptr(vert_x, 'XV', this%input_mempath)
    call mem_setptr(vert_y, 'YV', this%input_mempath)
    !
    ! -- set vertices 2d array
    if (associated(vert_x) .and. associated(vert_y)) then
      do i = 1, this%nvert
        this%vertices(1, i) = vert_x(i)
        this%vertices(2, i) = vert_y(i)
      end do
    else
      call store_error('Required Vertex arrays not found.')
    end if
    !
    ! -- log
    if (this%iout > 0) then
      write (this%iout, '(1x,a)') 'Discretization Vertex data loaded'
    end if
    !
    call memorystore_release('XV', this%input_mempath)
    call memorystore_release('YV', this%input_mempath)
  end subroutine source_vertices
 
  !> @brief Build data structures to hold cell vertex info
  !<
  subroutine define_cellverts(this, icell2d, ncvert, icvert)
    ! -- modules
    use sparsemodule, only: sparsematrix
    ! -- dummy
    class(disutype) :: this
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icell2d
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: ncvert
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icvert
    ! -- locals
    type(sparsematrix) :: vert_spm
    integer(I4B) :: i, j, ierr
    integer(I4B) :: icv_idx, startvert, maxnnz = 5
    !
    ! -- initialize sparse matrix
    call vert_spm%init(this%nodesuser, this%nvert, maxnnz)
    !
    ! -- add sparse matrix connections from input memory paths
    icv_idx = 1
    do i = 1, this%nodesuser
      if (icell2d(i) /= i) call store_error('ICELL2D input sequence violation.')
      do j = 1, ncvert(i)
        call vert_spm%addconnection(i, icvert(icv_idx), 0)
        if (j == 1) then
          startvert = icvert(icv_idx)
        elseif (j == ncvert(i) .and. (icvert(icv_idx) /= startvert)) then
          call vert_spm%addconnection(i, startvert, 0)
        end if
        icv_idx = icv_idx + 1
      end do
    end do
    !
    ! -- allocate and fill iavert and javert
    call mem_allocate(this%iavert, this%nodesuser + 1, 'IAVERT', this%memoryPath)
    call mem_allocate(this%javert, vert_spm%nnz, 'JAVERT', this%memoryPath)
    call vert_spm%filliaja(this%iavert, this%javert, ierr)
    call vert_spm%destroy()
    !
  end subroutine define_cellverts
 
  !> @brief Copy cell2d data from IDM into package
  !<
  subroutine source_cell2d(this)
    ! -- dummy
    class(disutype) :: this
    ! -- locals
    integer(I4B), dimension(:), contiguous, pointer :: icell2d => null()
    integer(I4B), dimension(:), contiguous, pointer :: ncvert => null()
    integer(I4B), dimension(:), contiguous, pointer :: icvert => null()
    real(DP), dimension(:), contiguous, pointer :: cell_x => null()
    real(DP), dimension(:), contiguous, pointer :: cell_y => null()
    integer(I4B) :: i
    !
    ! -- set pointers to input path ncvert and icvert
    call mem_setptr(icell2d, 'ICELL2D', this%input_mempath)
    call mem_setptr(ncvert, 'NCVERT', this%input_mempath)
    call mem_setptr(icvert, 'ICVERT', this%input_mempath)
    !
    ! --
    if (associated(icell2d) .and. associated(ncvert) &
        .and. associated(icvert)) then
      call this%define_cellverts(icell2d, ncvert, icvert)
    else
      call store_error('Required cell vertex arrays not found.')
    end if
    !
    ! -- set pointers to cell center arrays
    call mem_setptr(cell_x, 'XC', this%input_mempath)
    call mem_setptr(cell_y, 'YC', this%input_mempath)
    !
    ! -- set cell centers
    if (associated(cell_x) .and. associated(cell_y)) then
      do i = 1, this%nodesuser
        this%cellxy(1, i) = cell_x(i)
        this%cellxy(2, i) = cell_y(i)
      end do
    else
      call store_error('Required cell center arrays not found.')
    end if
    !
    ! -- log
    if (this%iout > 0) then
      write (this%iout, '(1x,a)') 'Discretization Cell2d data loaded'
    end if
    !
    call memorystore_release('ICELL2D', this%input_mempath)
    call memorystore_release('NCVERT', this%input_mempath)
    call memorystore_release('ICVERT', this%input_mempath)
    call memorystore_release('XC', this%input_mempath)
    call memorystore_release('YC', this%input_mempath)
  end subroutine source_cell2d
 
  !> @brief Write a binary grid file
  !<
  subroutine write_grb(this, icelltype)
    ! -- modules
    use openspecmodule, only: access, form
    use constantsmodule, only: lenbigline
    ! -- dummy
    class(disutype) :: this
    integer(I4B), dimension(:), intent(in) :: icelltype
    ! -- local
    integer(I4B) :: i, iunit, ntxt, version
    integer(I4B), parameter :: lentxt = 100
    character(len=50) :: txthdr
    character(len=lentxt) :: txt
    character(len=LINELENGTH) :: fname
    character(len=LENBIGLINE) :: crs
    logical(LGP) :: found_crs
    ! -- formats
    character(len=*), parameter :: fmtgrdsave = &
      "(4X,'BINARY GRID INFORMATION WILL BE WRITTEN TO:', &
       &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)"
    !
    ! -- Initialize
    version = 1
    ntxt = 11
    if (this%nvert > 0) ntxt = ntxt + 5
    !
    call mem_set_value(crs, 'CRS', this%input_mempath, found_crs)
    !
    ! -- set version
    if (found_crs) then
      ntxt = ntxt + 1
      version = 2
    end if
    !
    ! -- Open the file
    fname = trim(this%output_fname)
    iunit = getunit()
    write (this%iout, fmtgrdsave) iunit, trim(adjustl(fname))
    call openfile(iunit, this%iout, trim(adjustl(fname)), 'DATA(BINARY)', &
                  form, access, 'REPLACE')
    !
    ! -- write header information
    write (txthdr, '(a)') 'GRID DISU'
    txthdr(50:50) = new_line('a')
    write (iunit) txthdr
    write (txthdr, '(a, i0)') 'VERSION ', version
    txthdr(50:50) = new_line('a')
    write (iunit) txthdr
    write (txthdr, '(a, i0)') 'NTXT ', ntxt
    txthdr(50:50) = new_line('a')
    write (iunit) txthdr
    write (txthdr, '(a, i0)') 'LENTXT ', lentxt
    txthdr(50:50) = new_line('a')
    write (iunit) txthdr
    !
    ! -- write variable definitions
    write (txt, '(3a, i0)') 'NODES ', 'INTEGER ', 'NDIM 0 # ', this%nodesuser
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, i0)') 'NJA ', 'INTEGER ', 'NDIM 0 # ', this%con%nja
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, 1pg24.15)') 'XORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%xorigin
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, 1pg24.15)') 'YORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%yorigin
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, 1pg24.15)') 'ANGROT ', 'DOUBLE ', 'NDIM 0 # ', this%angrot
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, i0)') 'TOP ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, i0)') 'BOT ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, i0)') 'IA ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + 1
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, i0)') 'JA ', 'INTEGER ', 'NDIM 1 ', this%con%nja
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, i0)') 'IDOMAIN ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    write (txt, '(3a, i0)') 'ICELLTYPE ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
    txt(lentxt:lentxt) = new_line('a')
    write (iunit) txt
    !
    ! -- if vertices have been read then write additional header information
    if (this%nvert > 0) then
      write (txt, '(3a, i0)') 'VERTICES ', 'DOUBLE ', 'NDIM 2 2 ', this%nvert
      txt(lentxt:lentxt) = new_line('a')
      write (iunit) txt
      write (txt, '(3a, i0)') 'CELLX ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser
      txt(lentxt:lentxt) = new_line('a')
      write (iunit) txt
      write (txt, '(3a, i0)') 'CELLY ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser
      txt(lentxt:lentxt) = new_line('a')
      write (iunit) txt
      write (txt, '(3a, i0)') 'IAVERT ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + 1
      txt(lentxt:lentxt) = new_line('a')
      write (iunit) txt
      write (txt, '(3a, i0)') 'JAVERT ', 'INTEGER ', 'NDIM 1 ', size(this%javert)
      txt(lentxt:lentxt) = new_line('a')
      write (iunit) txt
    end if
    !
    ! -- if version 2 write character array headers
    if (version == 2) then
      if (found_crs) then
        write (txt, '(3a, i0)') 'CRS ', 'CHARACTER ', 'NDIM 1 ', &
          len_trim(crs)
        txt(lentxt:lentxt) = new_line('a')
        write (iunit) txt
      end if
    end if
    !
    ! -- write data
    write (iunit) this%nodesuser ! nodes
    write (iunit) this%nja ! nja
    write (iunit) this%xorigin ! xorigin
    write (iunit) this%yorigin ! yorigin
    write (iunit) this%angrot ! angrot
    write (iunit) this%top1d ! top
    write (iunit) this%bot1d ! bot
    write (iunit) this%con%iausr ! ia
    write (iunit) this%con%jausr ! ja
    write (iunit) this%idomain ! idomain
    write (iunit) icelltype ! icelltype
    !
    ! -- if vertices have been read then write additional data
    if (this%nvert > 0) then
      write (iunit) this%vertices ! vertices
      write (iunit) (this%cellxy(1, i), i=1, this%nodesuser) ! cellx
      write (iunit) (this%cellxy(2, i), i=1, this%nodesuser) ! celly
      write (iunit) this%iavert ! iavert
      write (iunit) this%javert ! javert
    end if
    !
    ! -- if version 2 write character array data
    if (version == 2) then
      if (found_crs) write (iunit) trim(crs) ! crs user input
    end if
    !
    ! -- Close the file
    close (iunit)
    !
  end subroutine write_grb
 
  !> @brief Get reduced node number from user node number
  !<
  function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber)
    class(disutype), intent(in) :: this
    integer(I4B), intent(in) :: nodeu
    integer(I4B), intent(in) :: icheck
    integer(I4B) :: nodenumber
    !
    if (icheck /= 0) then
      if (nodeu < 1 .or. nodeu > this%nodesuser) then
        write (errmsg, '(a,i0,a,i0,a)') &
          'Node number (', nodeu, ') is less than 1 or greater than nodes (', &
          this%nodesuser, ').'
        call store_error(errmsg)
      end if
    end if
    !
    ! -- set node number to passed in nodenumber since there is a one to one
    !    mapping for an unstructured grid
    if (this%nodes == this%nodesuser) then
      nodenumber = nodeu
    else
      nodenumber = this%nodereduced(nodeu)
    end if
    !
  end function get_nodenumber_idx1
 
  !> @brief Get normal vector components between the cell and a given neighbor
  !!
  !! The normal points outward from the shared face between noden and nodem.
  !<
  subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, &
                               ipos)
    ! -- dummy
    class(disutype) :: this
    integer(I4B), intent(in) :: noden !< cell (reduced nn)
    integer(I4B), intent(in) :: nodem !< neighbor (reduced nn)
    integer(I4B), intent(in) :: ihc !< horizontal connection flag
    real(DP), intent(inout) :: xcomp
    real(DP), intent(inout) :: ycomp
    real(DP), intent(inout) :: zcomp
    integer(I4B), intent(in) :: ipos
    ! -- local
    real(DP) :: angle, dmult
    !
    ! -- Set vector components based on ihc
    if (ihc == 0) then
      !
      ! -- connection is vertical
      xcomp = dzero
      ycomp = dzero
      if (nodem < noden) then
        !
        ! -- nodem must be above noden, so upward connection
        zcomp = done
      else
        !
        ! -- nodem must be below noden, so downward connection
        zcomp = -done
      end if
    else
      ! -- find from anglex, since anglex is symmetric, need to flip vector
      !    for lower triangle (nodem < noden)
      angle = this%con%anglex(this%con%jas(ipos))
      dmult = done
      if (nodem < noden) dmult = -done
      xcomp = cos(angle) * dmult
      ycomp = sin(angle) * dmult
      zcomp = dzero
    end if
    !
  end subroutine connection_normal
 
  !> @brief Get unit vector components between the cell and a given neighbor
  !!
  !! Saturation must be provided to compute cell center vertical coordinates.
  !! Also return the straight-line connection length.
  !<
  subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, &
                               xcomp, ycomp, zcomp, conlen)
    ! -- dummy
    class(disutype) :: this
    integer(I4B), intent(in) :: noden
    integer(I4B), intent(in) :: nodem
    logical, intent(in) :: nozee
    real(DP), intent(in) :: satn
    real(DP), intent(in) :: satm
    integer(I4B), intent(in) :: ihc
    real(DP), intent(inout) :: xcomp
    real(DP), intent(inout) :: ycomp
    real(DP), intent(inout) :: zcomp
    real(DP), intent(inout) :: conlen
    ! -- local
    real(DP) :: xn, xm, yn, ym, zn, zm
    !
    ! -- Terminate with error if requesting unit vector components for problems
    !    without cell data
    if (size(this%cellxy, 2) < 1) then
      write (errmsg, '(a)') &
        'Cannot calculate unit vector components for DISU grid if VERTEX '// &
        'data are not specified'
      call store_error(errmsg, terminate=.true.)
    end if
    !
    ! -- get xy center coords
    xn = this%xc(noden)
    yn = this%yc(noden)
    xm = this%xc(nodem)
    ym = this%yc(nodem)
    !
    ! -- Set vector components based on ihc
    if (ihc == 0) then
      !
      ! -- vertical connection, calculate z as cell center elevation
      zn = this%bot(noden) + dhalf * (this%top(noden) - this%bot(noden))
      zm = this%bot(nodem) + dhalf * (this%top(nodem) - this%bot(nodem))
    else
      !
      ! -- horizontal connection, with possible z component due to cell offsets
      !    and/or water table conditions
      if (nozee) then
        zn = dzero
        zm = dzero
      else
        zn = this%bot(noden) + dhalf * satn * (this%top(noden) - this%bot(noden))
        zm = this%bot(nodem) + dhalf * satm * (this%top(nodem) - this%bot(nodem))
      end if
    end if
    !
    ! -- Use coords to find vector components and connection length
    call line_unit_vector(xn, yn, zn, xm, ym, zm, xcomp, ycomp, zcomp, &
                          conlen)
    !
  end subroutine connection_vector
 
  !> @brief Get the discretization type
  !<
  subroutine get_dis_type(this, dis_type)
    ! -- dummy
    class(disutype), intent(in) :: this
    character(len=*), intent(out) :: dis_type
    !
    dis_type = "DISU"
    !
  end subroutine get_dis_type
 
  !> @brief Get the discretization type enumeration
  function get_dis_enum(this) result(dis_enum)
    use constantsmodule, only: disu
    class(disutype), intent(in) :: this
    integer(I4B) :: dis_enum
    dis_enum = disu
  end function get_dis_enum
 
  !> @brief Allocate and initialize scalar variables
  !<
  subroutine allocate_scalars(this, name_model, input_mempath)
    ! -- dummy
    class(disutype) :: this
    character(len=*), intent(in) :: name_model
    character(len=*), intent(in) :: input_mempath
    !
    ! -- Allocate parent scalars
    call this%DisBaseType%allocate_scalars(name_model, input_mempath)
    !
    ! -- Allocate variables for DISU
    call mem_allocate(this%njausr, 'NJAUSR', this%memoryPath)
    call mem_allocate(this%nvert, 'NVERT', this%memoryPath)
    call mem_allocate(this%nangldegxerr, 'NANGLDEGXERR', this%memoryPath)
    call mem_allocate(this%voffsettol, 'VOFFSETTOL', this%memoryPath)
    call mem_allocate(this%iangledegx, 'IANGLEDEGX', this%memoryPath)
    !
    ! -- Set values
    this%ndim = 1
    this%njausr = 0
    this%nvert = 0
    this%nangldegxerr = 0
    this%voffsettol = dzero
    this%iangledegx = 0
    this%readFromFile = .false.
    !
  end subroutine allocate_scalars
 
  !> @brief Allocate and initialize arrays
  !<
  subroutine allocate_arrays(this)
    ! -- dummy
    class(disutype) :: this
    !
    ! -- Allocate arrays in DisBaseType (mshape, top, bot, area)
    call this%DisBaseType%allocate_arrays()
    !
    ! -- Allocate arrays in DISU
    if (this%nodes < this%nodesuser) then
      call mem_allocate(this%nodeuser, this%nodes, 'NODEUSER', this%memoryPath)
      call mem_allocate(this%nodereduced, this%nodesuser, 'NODEREDUCED', &
                        this%memoryPath)
    else
      call mem_allocate(this%nodeuser, 1, 'NODEUSER', this%memoryPath)
      call mem_allocate(this%nodereduced, 1, 'NODEREDUCED', this%memoryPath)
    end if
    !
    ! -- Initialize
    this%mshape(1) = this%nodesuser
    !
  end subroutine allocate_arrays
 
  !> @brief Allocate arrays in memory manager
  !<
  subroutine allocate_arrays_mem(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(disutype) :: this
    !
    call mem_allocate(this%idomain, this%nodes, 'IDOMAIN', this%memoryPath)
    call mem_allocate(this%vertices, 2, this%nvert, 'VERTICES', this%memoryPath)
    if (this%icondir > 0) then
      call mem_allocate(this%cellxy, 2, this%nodes, 'CELLXY', this%memoryPath)
    else
      call mem_allocate(this%cellxy, 2, 0, 'CELLXY', this%memoryPath)
    end if
    !
  end subroutine allocate_arrays_mem
 
  !> @brief Convert a string to a user nodenumber
  !!
  !! Parse and return user nodenumber.
  !! If flag_string is present and true, the first token may be
  !! non-numeric (e.g. boundary name). In this case, return -2.
  !<
  function nodeu_from_string(this, lloc, istart, istop, in, iout, line, &
                             flag_string, allow_zero) result(nodeu)
    ! -- dummy
    class(disutype) :: this
    integer(I4B), intent(inout) :: lloc
    integer(I4B), intent(inout) :: istart
    integer(I4B), intent(inout) :: istop
    integer(I4B), intent(in) :: in
    integer(I4B), intent(in) :: iout
    character(len=*), intent(inout) :: line
    logical, optional, intent(in) :: flag_string
    logical, optional, intent(in) :: allow_zero
    integer(I4B) :: nodeu
    ! -- local
    integer(I4B) :: lloclocal, ndum, istat, n
    real(dp) :: r
    !
    if (present(flag_string)) then
      if (flag_string) then
        ! Check to see if first token in line can be read as an integer.
        lloclocal = lloc
        call urword(line, lloclocal, istart, istop, 1, ndum, r, iout, in)
        read (line(istart:istop), *, iostat=istat) n
        if (istat /= 0) then
          ! First token in line is not an integer; return flag to this effect.
          nodeu = -2
          return
        end if
      end if
    end if
    !
    call urword(line, lloc, istart, istop, 2, nodeu, r, iout, in)
    !
    if (nodeu == 0) then
      if (present(allow_zero)) then
        if (allow_zero) then
          return
        end if
      end if
    end if
    !
    if (nodeu < 1 .or. nodeu > this%nodesuser) then
      write (errmsg, '(a,i0,a)') &
        "Node number in list (", nodeu, ") is outside of the grid. "// &
        "Cell number cannot be determined in line '"// &
        trim(adjustl(line))//"'."
      call store_error(errmsg)
      call store_error_unit(in)
    end if
    !
  end function nodeu_from_string
 
  !> @brief Convert a cellid string to a user nodenumber
  !!
  !! If flag_string is present and true, the first token may be
  !! non-numeric (e.g. boundary name). In this case, return -2.
  !!
  !! If allow_zero is present and true, and all indices are zero, the
  !! result can be zero. If allow_zero is false, a zero in any index is an error.
  !<
  function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, &
                             allow_zero) result(nodeu)
    ! -- return
    integer(I4B) :: nodeu
    ! -- dummy
    class(disutype) :: this
    character(len=*), intent(inout) :: cellid
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    logical, optional, intent(in) :: flag_string
    logical, optional, intent(in) :: allow_zero
    ! -- local
    integer(I4B) :: lloclocal, istart, istop, ndum, n
    integer(I4B) :: istat
    real(dp) :: r
    !
    if (present(flag_string)) then
      if (flag_string) then
        ! Check to see if first token in cellid can be read as an integer.
        lloclocal = 1
        call urword(cellid, lloclocal, istart, istop, 1, ndum, r, iout, inunit)
        read (cellid(istart:istop), *, iostat=istat) n
        if (istat /= 0) then
          ! First token in cellid is not an integer; return flag to this effect.
          nodeu = -2
          return
        end if
      end if
    end if
    !
    lloclocal = 1
    call urword(cellid, lloclocal, istart, istop, 2, nodeu, r, iout, inunit)
    !
    if (nodeu == 0) then
      if (present(allow_zero)) then
        if (allow_zero) then
          return
        end if
      end if
    end if
    !
    if (nodeu < 1 .or. nodeu > this%nodesuser) then
      write (errmsg, '(a,i0,a)') &
        "Cell number cannot be determined for cellid ("// &
        trim(adjustl(cellid))//") and results in a user "// &
        "node number (", nodeu, ") that is outside of the grid."
      call store_error(errmsg)
      call store_error_unit(inunit)
    end if
    !
  end function nodeu_from_cellid
 
  !> @brief Indicates whether the grid discretization supports layers
  !<
  logical function supports_layers(this)
    ! -- dummy
    class(disutype) :: this
    !
    supports_layers = .false.
    !
  end function supports_layers
 
  !> @brief Get number of cells per layer (total nodes since DISU isn't layered)
  !<
  function get_ncpl(this)
    ! -- return
    integer(I4B) :: get_ncpl
    ! -- dummy
    class(disutype) :: this
    !
    get_ncpl = this%nodesuser
    !
  end function get_ncpl
 
  !> @brief Read an integer array
  !<
  subroutine read_int_array(this, line, lloc, istart, istop, iout, in, &
                            iarray, aname)
    ! -- dummy
    class(disutype), intent(inout) :: this
    character(len=*), intent(inout) :: line
    integer(I4B), intent(inout) :: lloc
    integer(I4B), intent(inout) :: istart
    integer(I4B), intent(inout) :: istop
    integer(I4B), intent(in) :: in
    integer(I4B), intent(in) :: iout
    integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: iarray
    character(len=*), intent(in) :: aname
    ! -- local
    integer(I4B) :: nval
    integer(I4B), dimension(:), pointer, contiguous :: itemp
    !
    ! -- Point the temporary pointer array, which is passed to the reading
    !    subroutine.  The temporary array will point to ibuff if it is a
    !    reduced structured system, or to iarray if it is an unstructured
    !    model.
    if (this%nodes < this%nodesuser) then
      nval = this%nodesuser
      itemp => this%ibuff
    else
      nval = this%nodes
      itemp => iarray
    end if
    !
    ! -- Read the array
    ! -- Read unstructured input
    call readarray(in, itemp, aname, this%ndim, nval, iout, 0)
    !
    ! -- If reduced model, then need to copy from itemp(=>ibuff) to iarray
    if (this%nodes < this%nodesuser) then
      call this%fill_grid_array(itemp, iarray)
    end if
    !
  end subroutine read_int_array
 
  !> @brief Read a double precision array
  !<
  subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, &
                            darray, aname)
    ! -- dummy
    class(disutype), intent(inout) :: this
    character(len=*), intent(inout) :: line
    integer(I4B), intent(inout) :: lloc
    integer(I4B), intent(inout) :: istart
    integer(I4B), intent(inout) :: istop
    integer(I4B), intent(in) :: in
    integer(I4B), intent(in) :: iout
    real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray
    character(len=*), intent(in) :: aname
    ! -- local
    integer(I4B) :: nval
    real(DP), dimension(:), pointer, contiguous :: dtemp
    !
    ! -- Point the temporary pointer array, which is passed to the reading
    !    subroutine.  The temporary array will point to dbuff if it is a
    !    reduced structured system, or to darray if it is an unstructured
    !    model.
    if (this%nodes < this%nodesuser) then
      nval = this%nodesuser
      dtemp => this%dbuff
    else
      nval = this%nodes
      dtemp => darray
    end if
    !
    ! -- Read the array
    call readarray(in, dtemp, aname, this%ndim, nval, iout, 0)
    !
    ! -- If reduced model, then need to copy from dtemp(=>dbuff) to darray
    if (this%nodes < this%nodesuser) then
      call this%fill_grid_array(dtemp, darray)
    end if
    !
  end subroutine read_dbl_array
 
  !> @brief Record a double precision array
  !!
  !! The array is written to a formatted or unformatted external file
  !! depending on the arguments.
  !<
  subroutine record_array(this, darray, iout, iprint, idataun, aname, &
                          cdatafmp, nvaluesp, nwidthp, editdesc, dinact)
    ! -- dummy
    class(disutype), intent(inout) :: this
    real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record
    integer(I4B), intent(in) :: iout !< ascii output unit number
    integer(I4B), intent(in) :: iprint !< whether to print the array
    integer(I4B), intent(in) :: idataun !< binary output unit number
    character(len=*), intent(in) :: aname !< text descriptor
    character(len=*), intent(in) :: cdatafmp ! write format
    integer(I4B), intent(in) :: nvaluesp !< values per line
    integer(I4B), intent(in) :: nwidthp !< number width
    character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E)
    real(DP), intent(in) :: dinact !< double precision value for cells excluded from model domain
    ! -- local
    integer(I4B) :: k, ifirst
    integer(I4B) :: nlay
    integer(I4B) :: nrow
    integer(I4B) :: ncol
    integer(I4B) :: nval
    integer(I4B) :: nodeu, noder
    integer(I4B) :: istart, istop
    real(DP), dimension(:), pointer, contiguous :: dtemp
    ! -- formats
    character(len=*), parameter :: fmthsv = &
      "(1X,/1X,a,' WILL BE SAVED ON UNIT ',I4, &
       &' AT END OF TIME STEP',I5,', STRESS PERIOD ',I4)"
    !
    ! -- set variables
    nlay = 1
    nrow = 1
    ncol = this%mshape(1)
    !
    ! -- If this is a reduced model, then copy the values from darray into
    !    dtemp.
    if (this%nodes < this%nodesuser) then
      nval = this%nodes
      dtemp => this%dbuff
      do nodeu = 1, this%nodesuser
        noder = this%get_nodenumber(nodeu, 0)
        if (noder <= 0) then
          dtemp(nodeu) = dinact
          cycle
        end if
        dtemp(nodeu) = darray(noder)
      end do
    else
      nval = this%nodes
      dtemp => darray
    end if
    !
    ! -- Print to iout if iprint /= 0
    if (iprint /= 0) then
      istart = 1
      do k = 1, nlay
        istop = istart + nrow * ncol - 1
        call ulaprufw(ncol, nrow, kstp, kper, k, iout, dtemp(istart:istop), &
                      aname, cdatafmp, nvaluesp, nwidthp, editdesc)
        istart = istop + 1
      end do
    end if
    !
    ! -- Save array to an external file.
    if (idataun > 0) then
      ! -- write to binary file by layer
      ifirst = 1
      istart = 1
      do k = 1, nlay
        istop = istart + nrow * ncol - 1
        if (ifirst == 1) write (iout, fmthsv) &
          trim(adjustl(aname)), idataun, &
          kstp, kper
        ifirst = 0
        call ulasav(dtemp(istart:istop), aname, kstp, kper, &
                    pertim, totim, ncol, nrow, k, idataun)
        istart = istop + 1
      end do
    elseif (idataun < 0) then
      !
      ! -- write entire array as one record
      call ubdsv1(kstp, kper, aname, -idataun, dtemp, ncol, nrow, nlay, &
                  iout, delt, pertim, totim)
    end if
    !
  end subroutine record_array
 
  !> @brief Record list header for imeth=6
  !<
  subroutine record_srcdst_list_header(this, text, textmodel, textpackage, &
                                       dstmodel, dstpackage, naux, auxtxt, &
                                       ibdchn, nlist, iout)
    ! -- dummy
    class(disutype) :: this
    character(len=16), intent(in) :: text
    character(len=16), intent(in) :: textmodel
    character(len=16), intent(in) :: textpackage
    character(len=16), intent(in) :: dstmodel
    character(len=16), intent(in) :: dstpackage
    integer(I4B), intent(in) :: naux
    character(len=16), dimension(:), intent(in) :: auxtxt
    integer(I4B), intent(in) :: ibdchn
    integer(I4B), intent(in) :: nlist
    integer(I4B), intent(in) :: iout
    ! -- local
    integer(I4B) :: nlay, nrow, ncol
    !
    nlay = 1
    nrow = 1
    ncol = this%mshape(1)
    !
    ! -- Use ubdsv06 to write list header
    call ubdsv06(kstp, kper, text, textmodel, textpackage, dstmodel, dstpackage, &
                 ibdchn, naux, auxtxt, ncol, nrow, nlay, &
                 nlist, iout, delt, pertim, totim)
    !
  end subroutine record_srcdst_list_header
 
  !> @brief Cast base to DISU
  !<
  !> @brief Get a 2D array of polygon vertices, listed in
  !!
  !! clockwise order beginning with the lower left corner.
  !! The array is empty if the optional cell vertices were
  !! not provided.
  !<
  subroutine get_polyverts(this, ic, polyverts, closed)
    ! -- dummy
    class(disutype), intent(inout) :: this
    integer(I4B), intent(in) :: ic !< cell number (reduced)
    real(DP), allocatable, intent(out) :: polyverts(:, :) !< polygon vertices (column-major indexing)
    logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
    ! -- local
    integer(I4B) :: icu, iavert, nverts, m, j
    logical(LGP) :: lclosed
 
    ! count vertices
    nverts = this%get_npolyverts(ic)
    if (nverts == 0) then
      allocate (polyverts(2, 0))
      return
    end if
 
    ! check closed option
    if (.not. (present(closed))) then
      lclosed = .false.
    else
      lclosed = closed
    end if
 
    ! allocate vertices array
    if (lclosed) then
      allocate (polyverts(2, nverts + 1))
    else
      allocate (polyverts(2, nverts))
    end if
 
    ! set vertices
    icu = this%get_nodeuser(ic)
    iavert = this%iavert(icu)
    do m = 1, nverts
      j = this%javert(iavert - 1 + m)
      polyverts(:, m) = (/this%vertices(1, j), this%vertices(2, j)/)
    end do
 
    ! close if enabled
    if (lclosed) &
      polyverts(:, nverts + 1) = polyverts(:, 1)
 
  end subroutine
 
  !> @brief Get the number of cell polygon vertices, 0 if none are defined.
  function get_npolyverts(this, ic, closed) result(npolyverts)
    class(disutype), intent(inout) :: this
    integer(I4B), intent(in) :: ic !< cell number (reduced)
    logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
    integer(I4B) :: npolyverts
    ! local
    integer(I4B) :: icu
 
    ! the vertices and cell2d blocks are optional
    if (this%nvert < 1) then
      npolyverts = 0
      return
    end if
 
    icu = this%get_nodeuser(ic)
    npolyverts = this%iavert(icu + 1) - this%iavert(icu) - 1
    if (present(closed)) then
      if (closed) npolyverts = npolyverts + 1
    end if
  end function get_npolyverts
 
  !> @brief Get the maximum number of cell polygon vertices.
  function get_max_npolyverts(this, closed) result(max_npolyverts)
    class(disutype), intent(inout) :: this
    logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
    integer(I4B) :: max_npolyverts
    integer(I4B) :: ic
 
    max_npolyverts = 0
    do ic = 1, this%nodes
      max_npolyverts = max(max_npolyverts, this%get_npolyverts(ic, closed))
    end do
  end function get_max_npolyverts
 
  function castasdisutype(dis) result(disu)
    ! -- dummy
    class(*), pointer :: dis !< base pointer to DISU object
    ! -- return
    class(disutype), pointer :: disu !< the resulting DISU pointer
    !
    disu => null()
    select type (dis)
    class is (disutype)
      disu => dis
    end select
    !
  end function castasdisutype
 
end module disumodule