BoundaryPackageExt.f90

Go to the documentation of this file.

!> @brief This module contains the extended boundary package
!!
!! This module contains the extended boundary type that itself
!! should be extended by model boundary packages that have been
!! updated to source static and dynamic input data from the
!! input context.
!!
!<
module bndextmodule
 
  use kindmodule, only: dp, lgp, i4b
  use constantsmodule, only: lenmempath, lenboundname, lenauxname, linelength
  use obsmodule, only: obs_cr
  use simvariablesmodule, only: errmsg
  use simmodule, only: store_error, count_errors, store_error_filename
  use bndmodule, only: bndtype
  use geomutilmodule, only: get_node, get_ijk
 
  implicit none
 
  private
  public :: bndexttype
 
  !> @ brief BndExtType
  !!
  !!  Generic extended boundary package type.  This derived type can be
  !!  overridden to define concrete boundary package types that source
  !!  all input from the input context.
  !<
  type, extends(bndtype) :: bndexttype
    ! -- characters
    ! -- scalars
    integer(I4B), pointer :: iper
    logical(LGP), pointer :: readarraygrid
    logical(LGP), pointer :: readasarrays
    ! -- arrays
    integer(I4B), dimension(:, :), pointer, contiguous :: cellid => null() !< input user cellid list
    integer(I4B), dimension(:), pointer, contiguous :: nodeulist => null() !< input user nodelist
  contains
    procedure :: bnd_df => bndext_df
    procedure :: bnd_rp => bndext_rp
    procedure :: bnd_da => bndext_da
    procedure :: allocate_scalars => bndext_allocate_scalars
    procedure :: allocate_arrays => bndext_allocate_arrays
    procedure :: source_options
    procedure :: source_dimensions
    procedure :: log_options
    procedure :: cellid_to_nlist
    procedure :: nodeu_to_nlist
    procedure :: layarr_to_nlist
    procedure :: default_nodelist
    procedure :: check_cellid
    procedure :: write_lstfile
    procedure :: bound_value
    procedure :: bnd_rp_log => bndext_rp_log
  end type bndexttype
 
  !> @ brief BndExtFoundType
  !!
  !!  This type is used to simplify the tracking of common parameters
  !!  that are sourced from the input context.
  !<
  type bndextfoundtype
    logical :: naux = .false.
    logical :: ipakcb = .false.
    logical :: iprpak = .false.
    logical :: iprflow = .false.
    logical :: boundnames = .false.
    logical :: auxmultname = .false.
    logical :: inewton = .false.
    logical :: auxiliary = .false.
    logical :: maxbound = .false.
  end type bndextfoundtype
 
contains
 
  !> @ brief Define boundary package options and dimensions
    !!
    !!  Define base boundary package options and dimensions for
    !!  a model boundary package.
    !!
  !<
  subroutine bndext_df(this, neq, dis)
    ! -- modules
    use basedismodule, only: disbasetype
    use timearrayseriesmanagermodule, only: timearrayseriesmanagertype, &
                                            tasmanager_cr
    use timeseriesmanagermodule, only: timeseriesmanagertype, tsmanager_cr
    ! -- dummy variables
    class(bndexttype), intent(inout) :: this !< BndExtType object
    integer(I4B), intent(inout) :: neq !< number of equations
    class(disbasetype), pointer :: dis !< discretization object
    !
    ! -- set pointer to dis object for the model
    this%dis => dis
    !
    ! -- Create time series managers
    ! -- Not in use by this type but BndType uses and deallocates
    call tsmanager_cr(this%TsManager, this%iout)
    call tasmanager_cr(this%TasManager, dis, this%name_model, this%iout)
    !
    ! -- create obs package
    call obs_cr(this%obs, this%inobspkg)
    !
    ! -- Write information to model list file
    write (this%iout, 1) trim(this%filtyp), trim(adjustl(this%text)), &
      this%input_mempath
1   format(1x, /1x, a, ' -- ', a, ' PACKAGE, VERSION 8, 2/22/2014', &
           ' INPUT READ FROM MEMPATH: ', a)
    !
    ! -- source options
    call this%source_options()
    !
    ! -- Define time series managers
    call this%tsmanager%tsmanager_df()
    call this%tasmanager%tasmanager_df()
    !
    ! -- source dimensions
    call this%source_dimensions()
    !
    ! -- update package moffset for packages that add rows
    if (this%npakeq > 0) then
      this%ioffset = neq - this%dis%nodes
    end if
    !
    ! -- update neq
    neq = neq + this%npakeq
    !
    ! -- Store information needed for observations
    if (this%bnd_obs_supported()) then
      call this%obs%obs_df(this%iout, this%packName, this%filtyp, this%dis)
      call this%bnd_df_obs()
    end if
    !
    ! -- Call define_listlabel to construct the list label that is written
    !    when PRINT_INPUT option is used.
    call this%define_listlabel()
  end subroutine bndext_df
 
  subroutine bndext_rp(this)
    ! -- modules
    use tdismodule, only: kper
    use memorymanagermodule, only: mem_deallocate, mem_reallocate
    use memorymanagerextmodule, only: mem_set_value
    ! -- dummy variables
    class(bndexttype), intent(inout) :: this !< BndExtType object
    ! -- local variables
    logical(LGP) :: found
    integer(I4B) :: n
    !
    if (this%iper /= kper) return
 
    if (.not. this%readasarrays) then
      ! -- copy nbound from input context
      call mem_set_value(this%nbound, 'NBOUND', this%input_mempath, &
                         found)
    end if
 
    if (this%readarraygrid) then
      call this%nodeu_to_nlist()
    else if (this%readasarrays) then
      call this%layarr_to_nlist()
    else
      call this%cellid_to_nlist()
      !
      ! -- update boundname string list
      if (this%inamedbound /= 0) then
        do n = 1, size(this%boundname_cst)
          this%boundname(n) = this%boundname_cst(n)
        end do
      end if
    end if
  end subroutine bndext_rp
 
  !> @brief Write the input list to the listing file if requested
  !!
  !! Called from model control files after bnd_rp(), which ensures
  !! bound_value() dispatches to the correct derived type.
  !<
  subroutine bndext_rp_log(this)
    ! -- dummy
    class(bndexttype), intent(inout) :: this
    !
    if (this%iprpak /= 0) then
      call this%write_lstfile()
    end if
  end subroutine bndext_rp_log
 
  !> @ brief Deallocate package memory
  !<
  subroutine bndext_da(this)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate, mem_setptr
    ! -- dummy variables
    class(bndexttype) :: this !< BndExtType object
    !
    ! -- deallocate checkin paths
    call mem_deallocate(this%cellid, 'CELLID', this%memoryPath)
    call mem_deallocate(this%nodeulist, 'NODEULIST', this%memoryPath)
    call mem_deallocate(this%boundname_cst, 'BOUNDNAME_IDM', this%memoryPath)
    call mem_deallocate(this%auxvar, 'AUXVAR_IDM', this%memoryPath)
    !
    ! -- reassign pointers for base class _da
    call mem_setptr(this%boundname_cst, 'BOUNDNAME_CST', this%memoryPath)
    call mem_setptr(this%auxvar, 'AUXVAR', this%memoryPath)
    !
    ! -- scalars
    deallocate (this%readarraygrid)
    deallocate (this%readasarrays)
    nullify (this%readarraygrid)
    nullify (this%readasarrays)
    nullify (this%iper)
    !
    ! -- deallocate
    call this%BndType%bnd_da()
  end subroutine bndext_da
 
  !> @ brief Allocate package scalars
    !!
    !!  Allocate and initialize base boundary package scalars. This method
    !!  only needs to be overridden if additional scalars are defined
    !!  for a specific package.
    !!
  !<
  subroutine bndext_allocate_scalars(this)
    ! -- modules
    use memorymanagermodule, only: mem_setptr
    ! -- dummy variables
    class(bndexttype) :: this !< BndExtType object
    ! -- local variables
    !
    ! -- allocate base BndType scalars
    call this%BndType%allocate_scalars()
    !
    ! -- set IPER pointer
    call mem_setptr(this%iper, 'IPER', this%input_mempath)
 
    ! -- allocate internal scalars
    allocate (this%readarraygrid)
    allocate (this%readasarrays)
 
    ! -- initialize internal scalars
    this%readarraygrid = .false.
    this%readasarrays = .false.
  end subroutine bndext_allocate_scalars
 
  !> @ brief Allocate package arrays
    !!
    !!  Allocate and initialize base boundary package arrays. This method
    !!  only needs to be overridden if additional arrays are defined
    !!  for a specific package.
    !!
  !<
  subroutine bndext_allocate_arrays(this, nodelist, auxvar)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate, mem_setptr, mem_checkin
    ! -- dummy variables
    class(bndexttype) :: this !< BndExtType object
    ! -- local variables
    integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist !< package nodelist
    real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar !< package aux variable array
    !
    ! -- allocate base BndType arrays
    call this%BndType%allocate_arrays(nodelist, auxvar)
    !
    ! -- set input context pointers
    call mem_setptr(this%cellid, 'CELLID', this%input_mempath)
    call mem_setptr(this%nodeulist, 'NODEULIST', this%input_mempath)
    call mem_setptr(this%boundname_cst, 'BOUNDNAME', this%input_mempath)
    !
    ! -- checkin input context pointers
    call mem_checkin(this%cellid, 'CELLID', this%memoryPath, &
                     'CELLID', this%input_mempath)
    call mem_checkin(this%nodeulist, 'NODEULIST', this%memoryPath, &
                     'NODEULIST', this%input_mempath)
    call mem_checkin(this%boundname_cst, lenboundname, 'BOUNDNAME_IDM', &
                     this%memoryPath, 'BOUNDNAME', this%input_mempath)
    !
    if (present(auxvar)) then
      ! no-op
    else
      ! -- set auxvar input context pointer
      call mem_setptr(this%auxvar, 'AUXVAR', this%input_mempath)
      !
      ! -- checkin auxvar input context pointer
      call mem_checkin(this%auxvar, 'AUXVAR_IDM', this%memoryPath, &
                       'AUXVAR', this%input_mempath)
    end if
  end subroutine bndext_allocate_arrays
 
  !> @ brief Source package options from input context
  !<
  subroutine source_options(this)
    ! -- modules
    use memorymanagermodule, only: mem_reallocate, mem_setptr !, get_isize
    use memorymanagerextmodule, only: mem_set_value
    use inputoutputmodule, only: getunit, openfile
    use characterstringmodule, only: characterstringtype
    use sourcecommonmodule, only: filein_fname
    ! -- dummy variables
    class(bndexttype), intent(inout) :: this !< BndExtType object
    ! -- local variables
    type(bndextfoundtype) :: found
    logical(LGP) :: found_readarr
    character(len=LENAUXNAME) :: sfacauxname
    integer(I4B) :: n
    !
    ! -- update defaults with idm sourced values
    call mem_set_value(this%naux, 'NAUX', this%input_mempath, found%naux)
    call mem_set_value(this%ipakcb, 'IPAKCB', this%input_mempath, found%ipakcb)
    call mem_set_value(this%iprpak, 'IPRPAK', this%input_mempath, found%iprpak)
    call mem_set_value(this%iprflow, 'IPRFLOW', this%input_mempath, found%iprflow)
    call mem_set_value(this%inamedbound, 'BOUNDNAMES', this%input_mempath, &
                       found%boundnames)
    call mem_set_value(sfacauxname, 'AUXMULTNAME', this%input_mempath, &
                       found%auxmultname)
    call mem_set_value(this%inewton, 'INEWTON', this%input_mempath, found%inewton)
    call mem_set_value(this%readarraygrid, 'READARRAYGRID', this%input_mempath, &
                       found_readarr)
    call mem_set_value(this%readasarrays, 'READASARRAYS', this%input_mempath, &
                       found_readarr)
    !
    ! -- log found options
    call this%log_options(found, sfacauxname)
    !
    ! -- reallocate aux arrays if aux variables provided
    if (found%naux .and. this%naux > 0) then
      call mem_reallocate(this%auxname, lenauxname, this%naux, &
                          'AUXNAME', this%memoryPath)
      call mem_reallocate(this%auxname_cst, lenauxname, this%naux, &
                          'AUXNAME_CST', this%memoryPath)
      call mem_set_value(this%auxname_cst, 'AUXILIARY', this%input_mempath, &
                         found%auxiliary)
      !
      do n = 1, this%naux
        this%auxname(n) = this%auxname_cst(n)
      end do
    end if
    !
    ! -- save flows option active
    if (found%ipakcb) this%ipakcb = -1
    !
    ! -- auxmultname provided
    if (found%auxmultname) this%iauxmultcol = -1
    !
    !
    ! -- enforce 0 or 1 OBS6_FILENAME entries in option block
    if (filein_fname(this%obs%inputFilename, 'OBS6_FILENAME', &
                     this%input_mempath, this%input_fname)) then
      this%obs%active = .true.
      this%obs%inUnitObs = getunit()
      call openfile(this%obs%inUnitObs, this%iout, this%obs%inputFilename, 'OBS')
    end if
    !
    ! -- no newton specified
    if (found%inewton) this%inewton = 0
    !
    ! -- AUXMULTNAME was specified, so find column of auxvar that will be multiplier
    if (this%iauxmultcol < 0) then
      !
      ! -- Error if no aux variable specified
      if (this%naux == 0) then
        write (errmsg, '(a,2(1x,a))') &
          'AUXMULTNAME was specified as', trim(adjustl(sfacauxname)), &
          'but no AUX variables specified.'
        call store_error(errmsg)
      end if
      !
      ! -- Assign mult column
      this%iauxmultcol = 0
      do n = 1, this%naux
        if (sfacauxname == this%auxname(n)) then
          this%iauxmultcol = n
          exit
        end if
      end do
      !
      ! -- Error if aux variable cannot be found
      if (this%iauxmultcol == 0) then
        write (errmsg, '(a,2(1x,a))') &
          'AUXMULTNAME was specified as', trim(adjustl(sfacauxname)), &
          'but no AUX variable found with this name.'
        call store_error(errmsg)
      end if
    end if
    !
    if (this%readasarrays) then
      if (.not. this%dis%supports_layers()) then
        errmsg = 'READASARRAYS option is not compatible with selected'// &
                 ' discretization type.'
        call store_error(errmsg)
        call store_error_filename(this%input_fname)
      end if
    end if
    !
    ! -- terminate if errors were detected
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
  end subroutine source_options
 
  !> @ brief Log package options
  !<
  subroutine log_options(this, found, sfacauxname)
    ! -- modules
    ! -- dummy variables
    class(bndexttype), intent(inout) :: this !< BndExtType object
    type(bndextfoundtype), intent(in) :: found
    character(len=*), intent(in) :: sfacauxname
    ! -- local variables
    ! -- format
    character(len=*), parameter :: fmtreadasarrays = &
      &"(4x, 'PACKAGE INPUT WILL BE READ AS LAYER ARRAYS.')"
    character(len=*), parameter :: fmtreadarraygrid = &
      &"(4x, 'PACKAGE INPUT WILL BE READ AS GRID ARRAYS.')"
    character(len=*), parameter :: fmtflow = &
      &"(4x, 'FLOWS WILL BE SAVED TO BUDGET FILE SPECIFIED IN OUTPUT CONTROL')"
    !
    ! -- log found options
    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text)) &
      //' BASE OPTIONS'
    !
    if (this%readasarrays) then
      write (this%iout, fmtreadasarrays)
    end if
    !
    if (this%readarraygrid) then
      write (this%iout, fmtreadarraygrid)
    end if
    !
    if (found%ipakcb) then
      write (this%iout, fmtflow)
    end if
    !
    if (found%iprpak) then
      write (this%iout, '(4x,a)') &
        'LISTS OF '//trim(adjustl(this%text))//' CELLS WILL BE PRINTED.'
    end if
    !
    if (found%iprflow) then
      write (this%iout, '(4x,a)') trim(adjustl(this%text))// &
        ' FLOWS WILL BE PRINTED TO LISTING FILE.'
    end if
    !
    if (found%boundnames) then
      write (this%iout, '(4x,a)') trim(adjustl(this%text))// &
        ' BOUNDARIES HAVE NAMES IN LAST COLUMN.'
    end if
    !
    if (found%auxmultname) then
      write (this%iout, '(4x,a,a)') &
        'AUXILIARY MULTIPLIER NAME: ', sfacauxname
    end if
    !
    if (found%inewton) then
      write (this%iout, '(4x,a)') &
        'NEWTON-RAPHSON method disabled for unconfined cells'
    end if
    !
    ! -- close logging block
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' BASE OPTIONS'
  end subroutine log_options
 
  !> @ brief Source package dimensions from input context
  !<
  subroutine source_dimensions(this)
    use memorymanagerextmodule, only: mem_set_value
    ! -- dummy variables
    class(bndexttype), intent(inout) :: this !< BndExtType object
    ! -- local variables
    type(bndextfoundtype) :: found
    !
    if (this%readasarrays) then
      this%maxbound = this%dis%get_ncpl()
    else
      ! -- open dimensions logging block
      write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
        ' BASE DIMENSIONS'
 
      ! -- update defaults with idm sourced values
      call mem_set_value(this%maxbound, 'MAXBOUND', this%input_mempath, &
                         found%maxbound)
 
      write (this%iout, '(4x,a,i7)') 'MAXBOUND = ', this%maxbound
 
      ! -- close logging block
      write (this%iout, '(1x,a)') &
        'END OF '//trim(adjustl(this%text))//' BASE DIMENSIONS'
    end if
    !
    ! -- verify dimensions were set
    if (this%maxbound <= 0) then
      write (errmsg, '(a)') 'MAXBOUND must be an integer greater than zero.'
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end if
  end subroutine source_dimensions
 
  !> @ brief Update package nodelist
    !!
    !! Convert period updated cellids to node numbers.
    !!
  !<
  subroutine cellid_to_nlist(this)
    ! -- modules
    use simvariablesmodule, only: errmsg
    ! -- dummy
    class(bndexttype) :: this !< BndExtType object
    ! -- local
    integer(I4B), dimension(:), pointer :: cellid
    integer(I4B) :: n, nodeu, noder
    character(len=LINELENGTH) :: nodestr
    !
    ! -- update nodelist
    do n = 1, this%nbound
      !
      ! -- set cellid
      cellid => this%cellid(:, n)
      !
      ! -- ensure cellid is valid, store an error otherwise
      call this%check_cellid(n, cellid, this%dis%mshape, this%dis%ndim)
      !
      ! -- Determine user node number
      if (this%dis%ndim == 1) then
        nodeu = cellid(1)
      elseif (this%dis%ndim == 2) then
        nodeu = get_node(cellid(1), 1, cellid(2), &
                         this%dis%mshape(1), 1, &
                         this%dis%mshape(2))
      else
        nodeu = get_node(cellid(1), cellid(2), cellid(3), &
                         this%dis%mshape(1), &
                         this%dis%mshape(2), &
                         this%dis%mshape(3))
      end if
      !
      ! -- update the nodelist
      if (this%dis%nodes < this%dis%nodesuser) then
        ! -- convert user to reduced node numbers
        noder = this%dis%get_nodenumber(nodeu, 0)
        if (noder <= 0) then
          call this%dis%nodeu_to_string(nodeu, nodestr)
          write (errmsg, *) &
            ' Cell is outside active grid domain: '// &
            trim(adjustl(nodestr))
          call store_error(errmsg)
        end if
        this%nodelist(n) = noder
      else
        this%nodelist(n) = nodeu
      end if
    end do
    !
    ! -- exit if errors were found
    if (count_errors() > 0) then
      write (errmsg, *) count_errors(), ' errors encountered.'
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end if
  end subroutine cellid_to_nlist
 
  !> @ brief Update package nodelist
  !!
  !! Convert period user nodes to reduced nodes
  !!
  !<
  subroutine nodeu_to_nlist(this)
    ! -- modules
    use memorymanagermodule, only: mem_setptr
    ! -- dummy
    class(bndexttype) :: this !< BndExtType object
    integer(I4B) :: n, noder, nodeuser, ninactive
 
    ninactive = 0
 
    ! -- Set the nodelist
    do n = 1, this%nbound
      nodeuser = this%nodeulist(n)
      noder = this%dis%get_nodenumber(nodeuser, 0)
      if (noder > 0) then
        this%nodelist(n) = noder
      else
        ninactive = ninactive + 1
      end if
    end do
 
    ! update nbound
    this%nbound = this%nbound - ninactive
  end subroutine nodeu_to_nlist
 
  !> @brief Update the nodelist based on layer number variable input
  !!
  !! This is a module scoped routine to check for I<filtyp>
  !! input. If array input was provided, INI<filtyp> and I<filtyp>
  !! will be allocated in the input context.  If the read
  !! state variable INI<filtyp> is set to 1 during this period
  !! update, I<filtyp> input was read and is used here to update
  !! the nodelist.
  !!
  !<
  subroutine layarr_to_nlist(this)
    ! -- modules
    use memorymanagermodule, only: mem_setptr
    use constantsmodule, only: lenvarname
    ! -- dummy
    class(bndexttype) :: this !< BndExtType object
    character(len=LENVARNAME) :: ilayname, inilayname
    character(len=24) :: aname = '     LAYER OR NODE INDEX'
    ! -- local
    integer(I4B), dimension(:), contiguous, &
      pointer :: ilay => null()
    integer(I4B), pointer :: inilay => null()
    !
    ! set ilay and read state variable names
    ilayname = 'I'//trim(this%filtyp)
    inilayname = 'INI'//trim(this%filtyp)
    !
    ! -- set pointer to input context read state variable
    call mem_setptr(inilay, inilayname, this%input_mempath)
    !
    ! -- check read state
    if (inilay == 1) then
      ! -- ilay variable was read this period
      !
      ! -- set pointer to input context layer index variable
      call mem_setptr(ilay, ilayname, this%input_mempath)
      !
      ! -- update nodelist
      call this%dis%nlarray_to_nodelist(ilay, this%nodelist, this%maxbound, &
                                        this%nbound, aname)
    end if
  end subroutine layarr_to_nlist
 
  !> @brief Assign default nodelist when READASARRAYS is specified.
  !!
  !! Equivalent to reading layer number array as CONSTANT 1
  !<
  subroutine default_nodelist(this)
    ! -- dummy
    class(bndexttype) :: this
    ! -- local
    integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nodeu, noder, ipos
    !
    if (this%readasarrays) then
      !
      ! -- set variables
      if (this%dis%ndim == 3) then
        nlay = this%dis%mshape(1)
        nrow = this%dis%mshape(2)
        ncol = this%dis%mshape(3)
      elseif (this%dis%ndim == 2) then
        nlay = this%dis%mshape(1)
        nrow = 1
        ncol = this%dis%mshape(2)
      end if
      !
      ! -- Populate nodelist
      ipos = 1
      il = 1
      do ir = 1, nrow
        do ic = 1, ncol
          nodeu = get_node(il, ir, ic, nlay, nrow, ncol)
          noder = this%dis%get_nodenumber(nodeu, 0)
          this%nodelist(ipos) = noder
          ipos = ipos + 1
        end do
      end do
      !
      ! -- Assign nbound
      this%nbound = ipos - 1
    end if
  end subroutine default_nodelist
 
  !> @ brief Check for valid cellid
  !<
  subroutine check_cellid(this, ii, cellid, mshape, ndim)
    ! -- modules
    use simvariablesmodule, only: errmsg
    ! -- dummy
    class(bndexttype) :: this !< BndExtType object
    ! -- local
    integer(I4B), intent(in) :: ii
    integer(I4B), dimension(:), intent(in) :: cellid !< cellid
    integer(I4B), dimension(:), intent(in) :: mshape !< model shape
    integer(I4B), intent(in) :: ndim !< size of mshape
    character(len=20) :: cellstr, mshstr
    character(len=*), parameter :: fmterr = &
      "('List entry ',i0,' contains cellid ',a,' but this cellid is invalid &
      &for model with shape ', a)"
    character(len=*), parameter :: fmtndim1 = &
                                   "('(',i0,')')"
    character(len=*), parameter :: fmtndim2 = &
                                   "('(',i0,',',i0,')')"
    character(len=*), parameter :: fmtndim3 = &
                                   "('(',i0,',',i0,',',i0,')')"
    select case (ndim)
    case (1)
      !
      if (cellid(1) < 1 .or. cellid(1) > mshape(1)) then
        write (cellstr, fmtndim1) cellid(1)
        write (mshstr, fmtndim1) mshape(1)
        write (errmsg, fmterr) ii, trim(adjustl(cellstr)), trim(adjustl(mshstr))
        call store_error(errmsg)
      end if
      !
    case (2)
      !
      if (cellid(1) < 1 .or. cellid(1) > mshape(1) .or. &
          cellid(2) < 1 .or. cellid(2) > mshape(2)) then
        write (cellstr, fmtndim2) cellid(1), cellid(2)
        write (mshstr, fmtndim2) mshape(1), mshape(2)
        write (errmsg, fmterr) ii, trim(adjustl(cellstr)), trim(adjustl(mshstr))
        call store_error(errmsg)
      end if
      !
    case (3)
      !
      if (cellid(1) < 1 .or. cellid(1) > mshape(1) .or. &
          cellid(2) < 1 .or. cellid(2) > mshape(2) .or. &
          cellid(3) < 1 .or. cellid(3) > mshape(3)) then
        write (cellstr, fmtndim3) cellid(1), cellid(2), cellid(3)
        write (mshstr, fmtndim3) mshape(1), mshape(2), mshape(3)
        write (errmsg, fmterr) ii, trim(adjustl(cellstr)), trim(adjustl(mshstr))
        call store_error(errmsg)
      end if
      !
    case default
    end select
  end subroutine check_cellid
 
  !> @ brief Log package stress period input
    !!
    !! Log period based input. This routine requires a package specific
    !! bound_value() routine to report accurate bound values.
    !!
  !<
  subroutine write_lstfile(this)
    ! -- modules
    use constantsmodule, only: linelength, lenboundname, &
                               tableft, tabcenter, dzero
    use inputoutputmodule, only: ulstlb
    use tablemodule, only: tabletype, table_cr
    ! -- dummy
    class(bndexttype) :: this !< BndExtType object
    ! -- local
    character(len=10) :: cpos
    character(len=LINELENGTH) :: tag
    character(len=LINELENGTH), allocatable, dimension(:) :: words
    integer(I4B) :: ntabrows
    integer(I4B) :: ntabcols
    integer(I4B) :: ipos
    integer(I4B) :: ii, jj, i, j, k, nod
    integer(I4B) :: ldim
    integer(I4B) :: naux
    type(tabletype), pointer :: inputtab => null()
    ! -- formats
    character(len=LINELENGTH) :: fmtlstbn
    !
    ! -- Determine sizes
    ldim = this%ncolbnd
    naux = size(this%auxvar, 1)
    !
    ! -- dimension table
    ntabrows = this%nbound
    !
    ! -- start building format statement to parse this%label, which
    !    contains the column headers (except for boundname and auxnames)
    ipos = index(this%listlabel, 'NO.')
    if (ipos /= 0) then
      write (cpos, '(i10)') ipos + 3
      fmtlstbn = '(a'//trim(adjustl(cpos))
    else
      fmtlstbn = '(a7'
    end if
    ! -- sequence number, layer, row, and column.
    if (size(this%dis%mshape) == 3) then
      ntabcols = 4
      fmtlstbn = trim(fmtlstbn)//',a7,a7,a7'
      !
      ! -- sequence number, layer, and cell2d.
    else if (size(this%dis%mshape) == 2) then
      ntabcols = 3
      fmtlstbn = trim(fmtlstbn)//',a7,a7'
      !
      ! -- sequence number and node.
    else
      ntabcols = 2
      fmtlstbn = trim(fmtlstbn)//',a7'
    end if
    !
    ! -- Add fields for non-optional real values
    ntabcols = ntabcols + ldim
    do i = 1, ldim
      fmtlstbn = trim(fmtlstbn)//',a16'
    end do
    !
    ! -- Add field for boundary name
    if (this%inamedbound == 1) then
      ntabcols = ntabcols + 1
      fmtlstbn = trim(fmtlstbn)//',a16'
    end if
    !
    ! -- Add fields for auxiliary variables
    ntabcols = ntabcols + naux
    do i = 1, naux
      fmtlstbn = trim(fmtlstbn)//',a16'
    end do
    fmtlstbn = trim(fmtlstbn)//')'
    !
    ! -- allocate words
    allocate (words(ntabcols))
    !
    ! -- parse this%listlabel into words
    read (this%listlabel, fmtlstbn) (words(i), i=1, ntabcols)
    !
    ! -- initialize the input table object
    call table_cr(inputtab, ' ', ' ')
    call inputtab%table_df(ntabrows, ntabcols, this%iout)
    !
    ! -- add the columns
    ipos = 1
    call inputtab%initialize_column(words(ipos), 10, alignment=tabcenter)
    !
    ! -- discretization
    do i = 1, size(this%dis%mshape)
      ipos = ipos + 1
      call inputtab%initialize_column(words(ipos), 7, alignment=tabcenter)
    end do
    !
    ! -- non-optional variables
    do i = 1, ldim
      ipos = ipos + 1
      call inputtab%initialize_column(words(ipos), 16, alignment=tabcenter)
    end do
    !
    ! -- boundname
    if (this%inamedbound == 1) then
      ipos = ipos + 1
      tag = 'BOUNDNAME'
      call inputtab%initialize_column(tag, lenboundname, alignment=tableft)
    end if
    !
    ! -- aux variables
    do i = 1, naux
      call inputtab%initialize_column(this%auxname(i), 16, alignment=tabcenter)
    end do
    !
    ! -- Write the table
    do ii = 1, this%nbound
      call inputtab%add_term(ii)
      !
      ! -- discretization
      if (size(this%dis%mshape) == 3) then
        nod = this%nodelist(ii)
        call get_ijk(nod, this%dis%mshape(2), this%dis%mshape(3), &
                     this%dis%mshape(1), i, j, k)
        call inputtab%add_term(k)
        call inputtab%add_term(i)
        call inputtab%add_term(j)
      else if (size(this%dis%mshape) == 2) then
        nod = this%nodelist(ii)
        call get_ijk(nod, 1, this%dis%mshape(2), this%dis%mshape(1), i, j, k)
        call inputtab%add_term(k)
        call inputtab%add_term(j)
      else
        nod = this%nodelist(ii)
        call inputtab%add_term(nod)
      end if
      !
      ! -- non-optional variables
      do jj = 1, ldim
        call inputtab%add_term(this%bound_value(jj, ii))
      end do
      !
      ! -- boundname
      if (this%inamedbound == 1) then
        call inputtab%add_term(this%boundname(ii))
      end if
      !
      ! -- aux variables
      do jj = 1, naux
        call inputtab%add_term(this%auxvar(jj, ii))
      end do
    end do
    !
    ! -- deallocate the local variables
    call inputtab%table_da()
    deallocate (inputtab)
    nullify (inputtab)
    deallocate (words)
  end subroutine write_lstfile
 
  !> @ brief Return a bound value
    !!
    !!  Return a bound value associated with an ncolbnd index
    !!  and row.  This function should be overridden in the
    !!  derived package class.
    !!
  !<
  function bound_value(this, col, row) result(bndval)
    ! -- modules
    use constantsmodule, only: dnodata
    ! -- dummy variables
    class(bndexttype), intent(inout) :: this !< BndExtType object
    integer(I4B), intent(in) :: col
    integer(I4B), intent(in) :: row
    ! -- result
    real(dp) :: bndval
    !
    ! -- override this return value by redefining this
    !    routine in the derived package.
    bndval = dnodata
  end function bound_value
 
end module bndextmodule