diff --git a/src/Model/GroundWaterFlow/gwf3dis8.f90 b/src/Model/GroundWaterFlow/gwf3dis8.f90
index d165c275a3f..e7b6d62abde 100644
--- a/src/Model/GroundWaterFlow/gwf3dis8.f90
+++ b/src/Model/GroundWaterFlow/gwf3dis8.f90
@@ -44,7 +44,7 @@ module GwfDisModule
     procedure :: nodeu_from_string
     procedure :: nodeu_from_cellid
     procedure :: supports_layers
-    procedure :: get_ncpl
+    procedure :: get_ncpl => get_ncpl
     procedure :: connection_vector
     procedure :: connection_normal
     ! -- private
@@ -56,8 +56,8 @@ module GwfDisModule
     procedure :: log_griddata
     procedure :: grid_finalize
     procedure :: write_grb
-    procedure :: allocate_scalars
-    procedure :: allocate_arrays
+    procedure :: allocate_scalars => allocate_scalars_dis
+    procedure :: allocate_arrays => allocate_arrays_dis
     !
     ! -- Read a node-sized model array (reduced or not)
     procedure :: read_int_array
@@ -159,7 +159,7 @@ subroutine dis3d_da(this)
     call memorylist_remove(this%name_model, 'DIS', idm_context)
     !
     ! -- DisBaseType deallocate
-    call this%DisBaseType%dis_da()
+    call this%dis_da_default()
     !
     ! -- Deallocate scalars
     call mem_deallocate(this%nlay)
@@ -866,7 +866,7 @@ function get_nodenumber_idx3(this, k, i, j, icheck) &
     return
   end function get_nodenumber_idx3
 
-  subroutine allocate_scalars(this, name_model, input_mempath)
+  subroutine allocate_scalars_dis(this, name_model, input_mempath)
 ! ******************************************************************************
 ! allocate_scalars -- Allocate and initialize scalars
 ! ******************************************************************************
@@ -881,7 +881,7 @@ subroutine allocate_scalars(this, name_model, input_mempath)
 ! ------------------------------------------------------------------------------
     !
     ! -- Allocate parent scalars
-    call this%DisBaseType%allocate_scalars(name_model, input_mempath)
+    call this%allocate_scalars_default(name_model, input_mempath)
     !
     ! -- Allocate
     call mem_allocate(this%nlay, 'NLAY', this%memoryPath)
@@ -896,9 +896,9 @@ subroutine allocate_scalars(this, name_model, input_mempath)
     !
     ! -- Return
     return
-  end subroutine allocate_scalars
+  end subroutine allocate_scalars_dis
 
-  subroutine allocate_arrays(this)
+  subroutine allocate_arrays_dis(this)
 ! ******************************************************************************
 ! allocate_arrays -- Allocate arrays
 ! ******************************************************************************
@@ -912,7 +912,7 @@ subroutine allocate_arrays(this)
 ! ------------------------------------------------------------------------------
     !
     ! -- Allocate arrays in DisBaseType (mshape, top, bot, area)
-    call this%DisBaseType%allocate_arrays()
+    call this%allocate_arrays_default()
     !
     ! -- Allocate arrays for GwfDisType
     if (this%nodes < this%nodesuser) then
@@ -931,7 +931,7 @@ subroutine allocate_arrays(this)
     !
     ! -- Return
     return
-  end subroutine allocate_arrays
+  end subroutine allocate_arrays_dis
 
   function nodeu_from_string(this, lloc, istart, istop, in, iout, line, &
                              flag_string, allow_zero) result(nodeu)
diff --git a/src/Model/GroundWaterFlow/gwf3disu8.f90 b/src/Model/GroundWaterFlow/gwf3disu8.f90
index d33bf55eb2a..edbe97cce8c 100644
--- a/src/Model/GroundWaterFlow/gwf3disu8.f90
+++ b/src/Model/GroundWaterFlow/gwf3disu8.f90
@@ -55,12 +55,12 @@ module GwfDisuModule
     procedure :: connection_normal
     procedure :: connection_vector
     procedure :: supports_layers
-    procedure :: get_ncpl
+    procedure :: get_ncpl => get_ncpl
     procedure, public :: record_array
     procedure, public :: record_srcdst_list_header
     ! -- private
-    procedure :: allocate_scalars
-    procedure :: allocate_arrays
+    procedure :: allocate_scalars => allocate_scalars_disu
+    procedure :: allocate_arrays => allocate_arrays_disu
     procedure :: allocate_arrays_mem
     procedure :: source_options
     procedure :: source_dimensions
@@ -78,6 +78,9 @@ module GwfDisuModule
     ! -- Read a node-sized model array (reduced or not)
     procedure :: read_int_array
     procedure :: read_dbl_array
+    !
+    procedure :: nlarray_to_nodelist
+    procedure :: read_layer_array
   end type GwfDisuType
 
 contains
@@ -482,7 +485,7 @@ subroutine disu_da(this)
     call mem_deallocate(this%nodereduced)
     !
     ! -- DisBaseType deallocate
-    call this%DisBaseType%dis_da()
+    call this%dis_da_default()
     !
     ! -- Return
     return
@@ -1305,7 +1308,7 @@ subroutine get_dis_type(this, dis_type)
 
   end subroutine get_dis_type
 
-  subroutine allocate_scalars(this, name_model, input_mempath)
+  subroutine allocate_scalars_disu(this, name_model, input_mempath)
 ! ******************************************************************************
 ! allocate_scalars -- Allocate and initialize scalar variables in this class
 ! ******************************************************************************
@@ -1322,7 +1325,7 @@ subroutine allocate_scalars(this, name_model, input_mempath)
 ! ------------------------------------------------------------------------------
     !
     ! -- Allocate parent scalars
-    call this%DisBaseType%allocate_scalars(name_model, input_mempath)
+    call this%allocate_scalars_default(name_model, input_mempath)
     !
     ! -- Allocate variables for DISU
     call mem_allocate(this%njausr, 'NJAUSR', this%memoryPath)
@@ -1340,9 +1343,9 @@ subroutine allocate_scalars(this, name_model, input_mempath)
     !
     ! -- Return
     return
-  end subroutine allocate_scalars
+  end subroutine allocate_scalars_disu
 
-  subroutine allocate_arrays(this)
+  subroutine allocate_arrays_disu(this)
 ! ******************************************************************************
 ! allocate_arrays -- Read discretization information from file
 ! ******************************************************************************
@@ -1357,7 +1360,7 @@ subroutine allocate_arrays(this)
 ! ------------------------------------------------------------------------------
     !
     ! -- Allocate arrays in DisBaseType (mshape, top, bot, area)
-    call this%DisBaseType%allocate_arrays()
+    call this%allocate_arrays_default()
     !
     ! -- Allocate arrays in DISU
     if (this%nodes < this%nodesuser) then
@@ -1374,7 +1377,7 @@ subroutine allocate_arrays(this)
     !
     ! -- Return
     return
-  end subroutine allocate_arrays
+  end subroutine allocate_arrays_disu
 
   subroutine allocate_arrays_mem(this)
     use MemoryManagerModule, only: mem_allocate
@@ -1809,4 +1812,52 @@ function CastAsDisuType(dis) result(disu)
 
   end function CastAsDisuType
 
+! todo: are the below needed??? can they be removed fromm DiscretizationBase?
+
+  subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname)
+    ! -- modules
+    use SimModule, only: store_error
+    use ConstantsModule, only: LINELENGTH
+    ! -- dummy
+    class(GwfDisuType) :: this
+    integer(I4B), intent(in) :: maxbnd
+    integer(I4B), dimension(:), pointer, contiguous :: darray
+    integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
+    integer(I4B), intent(inout) :: nbound
+    character(len=*), intent(in) :: aname
+    !
+    errmsg = 'Programmer error: nlarray_to_nodelist called for DISU grid.'
+    call store_error(errmsg, terminate=.TRUE.)
+
+  end subroutine nlarray_to_nodelist
+
+  subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, &
+    icolbnd, aname, inunit, iout)
+  ! ******************************************************************************
+  ! read_layer_array -- Read a 2d double array into col icolbnd of darray.
+  !                     For cells that are outside of the active domain,
+  !                     do not copy the array value into darray.
+  ! ******************************************************************************
+  !
+  !    SPECIFICATIONS:
+  ! ------------------------------------------------------------------------------
+  ! -- modules
+  ! -- dummy
+  class(GwfDisuType) :: this
+  integer(I4B), intent(in) :: ncolbnd
+  integer(I4B), intent(in) :: maxbnd
+  integer(I4B), dimension(maxbnd) :: nodelist
+  real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray
+  integer(I4B), intent(in) :: icolbnd
+  character(len=*), intent(in) :: aname
+  integer(I4B), intent(in) :: inunit
+  integer(I4B), intent(in) :: iout
+  !
+  !
+  errmsg = 'Programmer error: read_layer_array called for DISU grid.'
+  call store_error(errmsg, terminate=.TRUE.)
+  !
+  ! -- return
+  end subroutine read_layer_array
+
 end module GwfDisuModule
diff --git a/src/Model/GroundWaterFlow/gwf3disv8.f90 b/src/Model/GroundWaterFlow/gwf3disv8.f90
index 505cd680e89..f9c4b6bb125 100644
--- a/src/Model/GroundWaterFlow/gwf3disv8.f90
+++ b/src/Model/GroundWaterFlow/gwf3disv8.f90
@@ -48,7 +48,7 @@ module GwfDisvModule
     procedure :: connection_normal
     procedure :: connection_vector
     procedure :: supports_layers
-    procedure :: get_ncpl
+    procedure :: get_ncpl => get_ncpl
     ! -- private
     procedure :: source_options
     procedure :: source_dimensions
@@ -62,8 +62,8 @@ module GwfDisvModule
     procedure :: grid_finalize
     procedure :: connect
     procedure :: write_grb
-    procedure :: allocate_scalars
-    procedure :: allocate_arrays
+    procedure :: allocate_scalars => allocate_scalars_disv
+    procedure :: allocate_arrays => allocate_arrays_disv
     procedure :: get_cell2d_area
     !
     procedure :: read_int_array
@@ -180,7 +180,7 @@ subroutine disv_da(this)
                            context=idm_context)
     !
     ! -- DisBaseType deallocate
-    call this%DisBaseType%dis_da()
+    call this%dis_da_default()
     !
     ! -- Deallocate scalars
     call mem_deallocate(this%nlay)
@@ -1234,7 +1234,7 @@ subroutine get_dis_type(this, dis_type)
 
   end subroutine get_dis_type
 
-  subroutine allocate_scalars(this, name_model, input_mempath)
+  subroutine allocate_scalars_disv(this, name_model, input_mempath)
 ! ******************************************************************************
 ! allocate_scalars -- Allocate and initialize scalars
 ! ******************************************************************************
@@ -1250,7 +1250,7 @@ subroutine allocate_scalars(this, name_model, input_mempath)
 ! ------------------------------------------------------------------------------
     !
     ! -- Allocate parent scalars
-    call this%DisBaseType%allocate_scalars(name_model, input_mempath)
+    call this%allocate_scalars_default(name_model, input_mempath)
     !
     ! -- Allocate
     call mem_allocate(this%nlay, 'NLAY', this%memoryPath)
@@ -1265,9 +1265,9 @@ subroutine allocate_scalars(this, name_model, input_mempath)
     !
     ! -- Return
     return
-  end subroutine allocate_scalars
+  end subroutine allocate_scalars_disv
 
-  subroutine allocate_arrays(this)
+  subroutine allocate_arrays_disv(this)
 ! ******************************************************************************
 ! allocate_arrays -- Allocate arrays
 ! ******************************************************************************
@@ -1281,7 +1281,7 @@ subroutine allocate_arrays(this)
 ! ------------------------------------------------------------------------------
     !
     ! -- Allocate arrays in DisBaseType (mshape, top, bot, area)
-    call this%DisBaseType%allocate_arrays()
+    call this%allocate_arrays_default()
     !
     ! -- Allocate arrays for GwfDisvType
     if (this%nodes < this%nodesuser) then
@@ -1298,7 +1298,7 @@ subroutine allocate_arrays(this)
     !
     ! -- Return
     return
-  end subroutine allocate_arrays
+  end subroutine allocate_arrays_disv
 
   function get_cell2d_area(this, icell2d) result(area)
 ! ******************************************************************************
diff --git a/src/Model/ModelUtilities/DiscretizationBase.f90 b/src/Model/ModelUtilities/DiscretizationBase.f90
index adcbe831379..3f170bcecfd 100644
--- a/src/Model/ModelUtilities/DiscretizationBase.f90
+++ b/src/Model/ModelUtilities/DiscretizationBase.f90
@@ -22,7 +22,9 @@ module BaseDisModule
   public :: DisBaseType
   public :: dis_transform_xy
 
-  type :: DisBaseType
+  ! Abstract base type for grid discretizations.
+  type, abstract :: DisBaseType
+
     character(len=LENMEMPATH) :: memoryPath !< path for memory allocation
     character(len=LENMEMPATH) :: input_mempath = '' !< input context mempath
     character(len=LENMODELNAME), pointer :: name_model => null() !< name of the model
@@ -53,11 +55,31 @@ module BaseDisModule
     integer(I4B), dimension(:), pointer, contiguous :: nodereduced => null() !< (size:nodesuser)contains reduced nodenumber (size 0 if not reduced); -1 means vertical pass through, 0 is idomain = 0
     integer(I4B), dimension(:), pointer, contiguous :: nodeuser => null() !< (size:nodes) given a reduced nodenumber, provide the user nodenumber (size 0 if not reduced)
   contains
-    procedure :: dis_df
-    procedure :: dis_ac
-    procedure :: dis_mc
-    procedure :: dis_ar
-    procedure :: dis_da
+
+    ! Fortran doesn't support virtual procedures in abstract types,
+    ! so we define '..._default' procedures here, to be called from
+    ! overriding procedures in child types, and bind them here (and
+    ! in children) to generically-named procedures for external use.
+    ! For discussion see:
+    ! https://fortran-lang.discourse.group/t/call-overridden-procedure-of-abstract-parent-type/590/15
+
+    ! lifecycle procedures
+    procedure(dis_df), deferred :: dis_df
+    procedure :: dis_ac_default
+    procedure :: dis_ac => dis_ac_default
+    procedure :: dis_mc_default
+    procedure :: dis_mc => dis_mc_default
+    procedure :: dis_ar_default
+    procedure :: dis_ar => dis_ar_default
+    procedure :: dis_da_default
+    procedure :: dis_da => dis_da_default
+
+    ! allocate scalar and array variables
+    procedure :: allocate_scalars_default
+    procedure :: allocate_scalars => allocate_scalars_default
+    procedure :: allocate_arrays_default
+    procedure :: allocate_arrays => allocate_arrays_default
+
     ! -- helper functions
     !
     ! -- get_nodenumber is an overloaded integer function that will always
@@ -71,74 +93,322 @@ module BaseDisModule
     procedure :: get_nodenumber_idx1
     procedure :: get_nodenumber_idx2
     procedure :: get_nodenumber_idx3
-    procedure :: get_nodeuser
-    procedure :: nodeu_to_string
-    procedure :: nodeu_to_array
-    procedure :: nodeu_from_string
-    procedure :: nodeu_from_cellid
+
+    ! conversion between user and reduced node numbers, and utilities
+    ! to parse user node numbers from lines of text and cellid strings
+    procedure :: get_nodeuser ! todo rename to noder_from_nodeu?
     procedure :: noder_from_string
     procedure :: noder_from_cellid
-    procedure :: connection_normal
-    procedure :: connection_vector
-    procedure :: get_dis_type
-    procedure :: supports_layers
-    procedure :: allocate_scalars
-    procedure :: allocate_arrays
-    procedure :: get_ncpl
+    procedure :: noder_to_string
+    procedure :: noder_to_array
+    procedure(nodeu_from_string), deferred :: nodeu_from_string
+    procedure(nodeu_from_cellid), deferred :: nodeu_from_cellid
+    procedure(nodeu_to_string), deferred :: nodeu_to_string
+    procedure(nodeu_to_array), deferred :: nodeu_to_array
+    procedure(nlarray_to_nodelist), deferred :: nlarray_to_nodelist
+
+    ! calculate normal and unit vector components for shared faces
+    procedure(connection_normal), deferred :: connection_normal
+    procedure(connection_vector), deferred :: connection_vector
+
+    ! characterization and properties of discretization
+    procedure(get_dis_type), deferred :: get_dis_type
+    procedure(supports_layers), deferred :: supports_layers
+
+    ! get cell-related information
+    procedure(get_ncpl), deferred :: get_ncpl
     procedure :: get_cell_volume
-    procedure :: write_grb
-    !
+    procedure :: get_area
+    procedure :: get_area_factor
+    procedure :: highest_active ! todo rename to get_highest_active?
+
+    ! write a binary grid file
+    procedure(write_grb), deferred :: write_grb
+
+    ! generic procedure to read an array from a line of text
+    generic :: read_grid_array => read_int_array, read_dbl_array
     procedure :: read_int_array
     procedure :: read_dbl_array
-    generic, public :: read_grid_array => read_int_array, read_dbl_array
-    procedure, public :: read_layer_array
+    
+    ! read a 2D double array from file
+    procedure(read_layer_array), deferred, public :: read_layer_array
+
+    ! generic procedure to fill arrays indexed by reduced node number
     procedure :: fill_int_array
     procedure :: fill_dbl_array
-    generic, public :: fill_grid_array => fill_int_array, fill_dbl_array
-    procedure, public :: read_list
-    !
-    procedure, public :: record_array
-    procedure, public :: record_connection_array
-    procedure, public :: noder_to_string
-    procedure, public :: noder_to_array
-    procedure, public :: record_srcdst_list_header
+    generic :: fill_grid_array => fill_int_array, fill_dbl_array
+
+    ! read from input files and write to output files
+    procedure :: read_list
+    procedure(record_srcdst_list_header), deferred :: record_srcdst_list_header
     procedure, private :: record_srcdst_list_entry
-    generic, public :: record_mf6_list_entry => record_srcdst_list_entry
-    procedure, public :: nlarray_to_nodelist
-    procedure, public :: highest_active
-    procedure, public :: get_area
-    procedure, public :: get_area_factor
+    procedure(record_array), deferred :: record_array
+    procedure :: record_connection_array
+    generic :: record_mf6_list_entry => record_srcdst_list_entry
 
   end type DisBaseType
 
-contains
+  abstract interface
+  !> @brief Get the discretization type.
+  subroutine get_dis_type(this, dis_type)
+    import DisBaseType
+    class(DisBaseType), intent(in) :: this
+    character(len=*), intent(out) :: dis_type
+  end subroutine get_dis_type
+  end interface
 
+  abstract interface
+  !> @brief Define the discretization.
   subroutine dis_df(this)
-! ******************************************************************************
-! dis_df -- Read discretization information from DISU input file
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- dummy
+    import DisBaseType
     class(DisBaseType) :: this
-! ------------------------------------------------------------------------------
-    !
-    call store_error('Program error: DisBaseType method dis_df not &
-                     &implemented.', terminate=.TRUE.)
-    !
-    ! -- Return
-    return
   end subroutine dis_df
+  end interface
 
-  subroutine dis_ac(this, moffset, sparse)
-! ******************************************************************************
-! dis_ac -- Add connections to sparse based on cell connectivity
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- modules
+  abstract interface
+  !> @brief Write binary grid file. Called from AR procedure.
+  subroutine write_grb(this, icelltype)
+    import DisBaseType
+    import I4B
+    class(DisBaseType) :: this
+    integer(I4B), dimension(:), intent(in) :: icelltype
+  end subroutine write_grb
+  end interface
+
+  abstract interface
+  !> @brief Convert user node number to a string in the form of (nodenumber) or (k,i,j).
+  subroutine nodeu_to_string(this, nodeu, str)
+    import DisBaseType
+    import I4B
+    class(DisBaseType) :: this
+    integer(I4B), intent(in) :: nodeu
+    character(len=*), intent(inout) :: str
+  end subroutine nodeu_to_string
+  end interface
+
+  abstract interface
+  !> @brief Convert user node number to  array with format (nodenumber) or (k,j) or (k,i,j).
+  subroutine nodeu_to_array(this, nodeu, arr)
+    import DisBaseType
+    import I4B
+    class(DisBaseType) :: this
+    integer(I4B), intent(in) :: nodeu
+    integer(I4B), dimension(:), intent(inout) :: arr
+  end subroutine nodeu_to_array
+  end interface
+
+  abstract interface
+  !> @brief Calculate normal vector components for shared cell faces.
+  !!
+  !! Reduced nodenumber cell is noden and its shared face with is cell nodem.
+  !! ihc is the horizontal connection flag.
+  subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, &
+                               ipos)
+    import DisBaseType
+    import I4B
+    import DP
+    class(DisBaseType) :: this
+    integer(I4B), intent(in) :: noden
+    integer(I4B), intent(in) :: nodem
+    integer(I4B), intent(in) :: ihc
+    real(DP), intent(inout) :: xcomp
+    real(DP), intent(inout) :: ycomp
+    real(DP), intent(inout) :: zcomp
+    integer(I4B), intent(in) :: ipos
+  end subroutine connection_normal
+  end interface
+
+  abstract interface
+  !> @brief Calculate unit vector components for shared cell faces.
+  !!
+  !! Reduced nodenumber cell is noden and its neighbor cell is nodem.
+  !! Saturations for these cells are also required so the vertical
+  !! position of the cell centers can be calculated. ihc is the
+  !! horizontal flag. Also return the straight-line connection length.
+  subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, &
+                               xcomp, ycomp, zcomp, conlen)
+    import DisBaseType
+    import I4B
+    import DP
+    class(DisBaseType) :: 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
+  end subroutine connection_vector
+  end interface
+
+  abstract interface
+  !> @brief Parse a user nodenumber from a line of text.
+  !!
+  !! The model is unstructured; just read user nodenumber.
+  !! If flag_string argument is present and true, the first token in string
+  !! is allowed to be a string (e.g. boundary name). In this case, if a string
+  !! is encountered, return value as -2.
+  function nodeu_from_string(this, lloc, istart, istop, in, iout, line, &
+                             flag_string, allow_zero) result(nodeu)
+    import DisBaseType
+    import I4B
+    class(DisBaseType) :: 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
+  end function nodeu_from_string
+  end interface
+
+  abstract interface
+  !> @brief Parse a user nodenumber from a cell ID string.
+  !!
+  !! If flag_string argument is present and true, the first token in string
+  !! is allowed to be a string (e.g. boundary name). In this case, if a string
+  !! is encountered, return value as -2.
+  !! If allow_zero argument is present and true, if all indices equal zero, the
+  !! result can be zero. If allow_zero is false, a zero in any index causes an
+  !! error.
+  function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, &
+                             allow_zero) result(nodeu)
+    import DisBaseType
+    import I4B
+    class(DisBaseType) :: 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
+    integer(I4B) :: nodeu
+  end function nodeu_from_cellid
+  end interface
+
+  abstract interface
+  !> @brief Indicates whether the discretization supports layers.
+  logical function supports_layers(this)
+    import DisBaseType
+    class(DisBaseType) :: this
+  end function supports_layers
+  end interface
+
+  abstract interface
+  !> @brief Get the number of cells per layer (nodes for DISU, since no layers).
+  function get_ncpl(this)
+    import DisBaseType
+    import I4B
+    integer(I4B) :: get_ncpl
+    class(DisBaseType) :: this
+  end function get_ncpl
+  end interface
+
+  abstract interface
+  !> @brief Read a 2d double array into col icolbnd of darray.
+  !! For cells that are outside of the active domain,
+  !! do not copy the array value into darray.
+  subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, &
+                              icolbnd, aname, inunit, iout)
+    import DisBaseType
+    import I4B
+    import DP
+    class(DisBaseType) :: this
+    integer(I4B), intent(in) :: ncolbnd
+    integer(I4B), intent(in) :: maxbnd
+    integer(I4B), dimension(maxbnd) :: nodelist
+    real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray
+    integer(I4B), intent(in) :: icolbnd
+    character(len=*), intent(in) :: aname
+    integer(I4B), intent(in) :: inunit
+    integer(I4B), intent(in) :: iout
+  end subroutine read_layer_array
+  end interface
+
+  abstract interface
+  !> @brief Record a double precision array.
+  !!
+  !! The array will be printed to an external file and/or written to an
+  !! unformatted external file depending on the argument specifications.
+  !!
+  !! SPECIFICATIONS:
+  !!   darray is the double precision array to record
+  !!   iout is the unit number for ascii output
+  !!   iprint is a flag indicating whether or not to print the array
+  !!   idataun is the unit number to which the array will be written in binary
+  !!     form; if negative then do not write by layers, write entire array
+  !!   aname is the text descriptor of the array
+  !!   cdatafmp is the fortran format for writing the array
+  !!   nvaluesp is the number of values per line for printing
+  !!   nwidthp is the width of the number for printing
+  !!   editdesc is the format type (I, G, F, S, E)
+  !!   dinact is the double precision value to use for cells that are excluded
+  !!     from the model domain
+  subroutine record_array(this, darray, iout, iprint, idataun, aname, &
+                          cdatafmp, nvaluesp, nwidthp, editdesc, dinact)
+    import DisBaseType
+    import I4B
+    import DP
+    class(DisBaseType), intent(inout) :: this
+    real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray
+    integer(I4B), intent(in) :: iout
+    integer(I4B), intent(in) :: iprint
+    integer(I4B), intent(in) :: idataun
+    character(len=*), intent(in) :: aname
+    character(len=*), intent(in) :: cdatafmp
+    integer(I4B), intent(in) :: nvaluesp
+    integer(I4B), intent(in) :: nwidthp
+    character(len=*), intent(in) :: editdesc
+    real(DP), intent(in) :: dinact
+  end subroutine record_array
+  end interface
+
+  abstract interface
+  !> @brief Record list header for imeth=6
+  subroutine record_srcdst_list_header(this, text, textmodel, textpackage, &
+                                       dstmodel, dstpackage, naux, auxtxt, &
+                                       ibdchn, nlist, iout)
+    import DisBaseType
+    import I4B
+    class(DisBaseType) :: 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
+  end subroutine record_srcdst_list_header
+  end interface
+
+  abstract interface
+  !> @brief Convert an integer array to nodelist.
+  !! For structured model, integer array is layer number; for unstructured
+  !! model, integer array is node number.
+  subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname)
+    import DisBaseType
+    import I4B
+    class(DisBaseType) :: this
+    integer(I4B), intent(in) :: maxbnd
+    integer(I4B), dimension(:), pointer, contiguous :: darray
+    integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
+    integer(I4B), intent(inout) :: nbound
+    character(len=*), intent(in) :: aname
+  end subroutine nlarray_to_nodelist
+  end interface
+
+contains
+
+  !> @brief Add connections to sparse based on cell connectivity.
+  subroutine dis_ac_default(this, moffset, sparse)
     use SparseModule, only: sparsematrix
     ! -- dummy
     class(DisBaseType) :: this
@@ -146,7 +416,6 @@ subroutine dis_ac(this, moffset, sparse)
     type(sparsematrix), intent(inout) :: sparse
     ! -- local
     integer(I4B) :: i, j, ipos, iglo, jglo
-! ------------------------------------------------------------------------------
     !
     do i = 1, this%nodes
       do ipos = this%con%ia(i), this%con%ia(i + 1) - 1
@@ -156,20 +425,10 @@ subroutine dis_ac(this, moffset, sparse)
         call sparse%addconnection(iglo, jglo, 1)
       end do
     end do
-    !
-    ! -- Return
-    return
-  end subroutine dis_ac
-
-  subroutine dis_mc(this, moffset, idxglo, matrix_sln)
-! ******************************************************************************
-! dis_mc -- Map the positions of cell connections in the numerical solution
-!   coefficient matrix.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- modules
+  end subroutine dis_ac_default
+
+  !> @brief Map cell connection positions in numerical solution coefficient matrix.
+  subroutine dis_mc_default(this, moffset, idxglo, matrix_sln)
     ! -- dummy
     class(DisBaseType) :: this
     integer(I4B), intent(in) :: moffset
@@ -177,7 +436,6 @@ subroutine dis_mc(this, moffset, idxglo, matrix_sln)
     class(MatrixBaseType), pointer :: matrix_sln
     ! -- local
     integer(I4B) :: i, j, ipos, iglo, jglo
-! ------------------------------------------------------------------------------
     !
     do i = 1, this%nodes
       iglo = i + moffset
@@ -187,26 +445,16 @@ subroutine dis_mc(this, moffset, idxglo, matrix_sln)
         idxglo(ipos) = matrix_sln%get_position(iglo, jglo)
       end do
     end do
-    !
-    ! -- Return
-    return
-  end subroutine dis_mc
-
-  subroutine dis_ar(this, icelltype)
-! ******************************************************************************
-! dis_ar -- Called from AR procedure.  Only task is to write binary grid file.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- modules
+  end subroutine dis_mc_default
+
+  !> @brief Allocate and read.
+  subroutine dis_ar_default(this, icelltype)
     ! -- dummy
     class(DisBaseType) :: this
     integer(I4B), dimension(:), intent(in) :: icelltype
     ! -- local
     integer(I4B), dimension(:), allocatable :: ict
     integer(I4B) :: nu, nr
-! ------------------------------------------------------------------------------
     !
     ! -- Expand icelltype to full grid; fill with 0 if cell is excluded
     allocate (ict(this%nodesuser))
@@ -220,45 +468,14 @@ subroutine dis_ar(this, icelltype)
     end do
     !
     if (this%nogrb == 0) call this%write_grb(ict)
-    !
-    ! -- Return
-    return
-  end subroutine dis_ar
-
-  subroutine write_grb(this, icelltype)
-! ******************************************************************************
-! write_grb -- Called from AR procedure.  Only task is to write binary grid file.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- modules
-    ! -- dummy
-    class(DisBaseType) :: this
-    integer(I4B), dimension(:), intent(in) :: icelltype
-    ! -- local
-! ------------------------------------------------------------------------------
-    !
-    !
-    call store_error('Program error: DisBaseType method write_grb not &
-                     &implemented.', terminate=.TRUE.)
-    !
-    ! -- Return
-    return
-  end subroutine write_grb
+  end subroutine dis_ar_default
 
-  subroutine dis_da(this)
-! ******************************************************************************
-! dis_da -- Deallocate discretization object
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
+  !> @brief Deallocate discretization object.
+  subroutine dis_da_default(this)
     ! -- modules
     use MemoryManagerModule, only: mem_deallocate
     ! -- dummy
     class(DisBaseType) :: this
-! ------------------------------------------------------------------------------
     !
     ! -- Strings
     deallocate (this%name_model)
@@ -292,68 +509,15 @@ subroutine dis_da(this)
     ! -- Connections
     call this%con%con_da()
     deallocate (this%con)
-    !
-    ! -- Return
-    return
-  end subroutine dis_da
-
-  subroutine nodeu_to_string(this, nodeu, str)
-! ******************************************************************************
-! nodeu_to_string -- Convert user node number to a string in the form of
-! (nodenumber) or (k,i,j)
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- dummy
-    class(DisBaseType) :: this
-    integer(I4B), intent(in) :: nodeu
-    character(len=*), intent(inout) :: str
-    ! -- local
-! ------------------------------------------------------------------------------
-    !
-    call store_error('Program error: DisBaseType method nodeu_to_string not &
-                     &implemented.', terminate=.TRUE.)
-    !
-    ! -- return
-    return
-  end subroutine nodeu_to_string
-
-  subroutine nodeu_to_array(this, nodeu, arr)
-! ******************************************************************************
-! nodeu_to_array -- Convert user node number to cellid and fill array with
-!                   (nodenumber) or (k,j) or (k,i,j)
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- dummy
-    class(DisBaseType) :: this
-    integer(I4B), intent(in) :: nodeu
-    integer(I4B), dimension(:), intent(inout) :: arr
-    ! -- local
-! ------------------------------------------------------------------------------
-    !
-    call store_error('Program error: DisBaseType method nodeu_to_array not &
-                     &implemented.', terminate=.TRUE.)
-    !
-    ! -- return
-    return
-  end subroutine nodeu_to_array
+  end subroutine dis_da_default
 
+  !> @brief Get user nodenumber from reduced node number.
   function get_nodeuser(this, noder) result(nodenumber)
-! ******************************************************************************
-! get_nodeuser -- Return the user nodenumber from the reduced node number
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- return
     integer(I4B) :: nodenumber
     ! -- dummy
     class(DisBaseType) :: this
     integer(I4B), intent(in) :: noder
-! ------------------------------------------------------------------------------
     !
     if (this%nodes < this%nodesuser) then
       nodenumber = this%nodeuser(noder)
@@ -365,16 +529,10 @@ function get_nodeuser(this, noder) result(nodenumber)
     return
   end function get_nodeuser
 
+  !> @brief Get nodenumber from the user specified node number.
+  !! Overridde in child classes to map to a model node number,
+  !! optionally performing a check.
   function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber)
-! ******************************************************************************
-! get_nodenumber -- Return a nodenumber from the user specified node number
-!                   with an option to perform a check.  This subroutine
-!                   can be overridden by child classes to perform mapping
-!                   to a model node number
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- modules
     use ConstantsModule, only: LINELENGTH
     use SimModule, only: store_error
@@ -384,7 +542,6 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber)
     integer(I4B), intent(in) :: icheck
     ! -- local
     integer(I4B) :: nodenumber
-! ------------------------------------------------------------------------------
     !
     nodenumber = 0
     call store_error('Program error: get_nodenumber_idx1 not implemented.', &
@@ -394,21 +551,14 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber)
     return
   end function get_nodenumber_idx1
 
+  !> @brief Override in child classes to map to a model node number.
   function get_nodenumber_idx2(this, k, j, icheck) result(nodenumber)
-! ******************************************************************************
-! get_nodenumber_idx2 -- This function should never be called.  It must be
-!   overridden by a child class.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     use SimModule, only: store_error
     ! -- dummy
     class(DisBaseType), intent(in) :: this
     integer(I4B), intent(in) :: k, j
     integer(I4B), intent(in) :: icheck
     integer(I4B) :: nodenumber
-! ------------------------------------------------------------------------------
     !
     nodenumber = 0
     call store_error('Program error: get_nodenumber_idx2 not implemented.', &
@@ -418,21 +568,15 @@ function get_nodenumber_idx2(this, k, j, icheck) result(nodenumber)
     return
   end function get_nodenumber_idx2
 
+  !> @brief This function will not be invoked for an unstructured
+  !! model, but it may be from a Discretization3dType model.
   function get_nodenumber_idx3(this, k, i, j, icheck) result(nodenumber)
-! ******************************************************************************
-! get_nodenumber_idx3 -- This function will not be invoked for an unstructured
-! model, but it may be from a Discretization3dType model.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     use SimModule, only: store_error
     ! -- dummy
     class(DisBaseType), intent(in) :: this
     integer(I4B), intent(in) :: k, i, j
     integer(I4B), intent(in) :: icheck
     integer(I4B) :: nodenumber
-! ------------------------------------------------------------------------------
     !
     nodenumber = 0
     call store_error('Program error: get_nodenumber_idx3 not implemented.', &
@@ -442,75 +586,7 @@ function get_nodenumber_idx3(this, k, i, j, icheck) result(nodenumber)
     return
   end function get_nodenumber_idx3
 
-  subroutine connection_normal(this, noden, nodem, ihc, xcomp, ycomp, zcomp, &
-                               ipos)
-! ******************************************************************************
-! connection_normal -- calculate the normal vector components for reduced
-!   nodenumber cell (noden) and its shared face with cell nodem.  ihc is the
-!   horizontal connection flag.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- modules
-    use SimModule, only: store_error
-    ! -- dummy
-    class(DisBaseType) :: this
-    integer(I4B), intent(in) :: noden
-    integer(I4B), intent(in) :: nodem
-    integer(I4B), intent(in) :: ihc
-    real(DP), intent(inout) :: xcomp
-    real(DP), intent(inout) :: ycomp
-    real(DP), intent(inout) :: zcomp
-    integer(I4B), intent(in) :: ipos
-! ------------------------------------------------------------------------------
-    !
-    call store_error('Program error: connection_normal not implemented.', &
-                     terminate=.TRUE.)
-    !
-    ! -- return
-    return
-  end subroutine connection_normal
-
-  subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, &
-                               xcomp, ycomp, zcomp, conlen)
-! ******************************************************************************
-! connection_vector -- calculate the unit vector components from reduced
-!   nodenumber cell (noden) to its neighbor cell (nodem).  The saturation for
-!   for these cells are also required so that the vertical position of the cell
-!   cell centers can be calculated.  ihc is the horizontal flag.  Also return
-!   the straight-line connection length.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- modules
-    use SimModule, only: store_error
-    ! -- dummy
-    class(DisBaseType) :: 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
-! ------------------------------------------------------------------------------
-    !
-    call store_error('Program error: connection_vector not implemented.', &
-                     terminate=.TRUE.)
-    !
-    ! -- return
-    return
-  end subroutine connection_vector
-
-  !> @brief get the x,y for a node transformed into
-  !! 'global coordinates' using xorigin, yorigin, angrot,
-  !< analogously to how flopy does this.
+  !> @brief Transform local to global coords using xorigin, yorigin, angrot.
   subroutine dis_transform_xy(x, y, xorigin, yorigin, angrot, xglo, yglo)
     real(DP), intent(in) :: x !< the cell-x coordinate to transform
     real(DP), intent(in) :: y !< the cell-y coordinate to transform
@@ -538,27 +614,8 @@ subroutine dis_transform_xy(x, y, xorigin, yorigin, angrot, xglo, yglo)
 
   end subroutine dis_transform_xy
 
-  !> @brief return discretization type
-  !<
-  subroutine get_dis_type(this, dis_type)
-    class(DisBaseType), intent(in) :: this
-    character(len=*), intent(out) :: dis_type
-
-    ! suppress warning
-    dis_type = "Not implemented"
-
-    call store_error('Program error: get_dis_type not implemented.', &
-                     terminate=.TRUE.)
-
-  end subroutine get_dis_type
-
-  subroutine allocate_scalars(this, name_model, input_mempath)
-! ******************************************************************************
-! allocate_scalars -- Allocate and initialize scalar variables in this class
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
+  !> @brief Allocate and initialize scalar variables.
+  subroutine allocate_scalars_default(this, name_model, input_mempath)
     ! -- modules
     use MemoryManagerModule, only: mem_allocate
     use MemoryManagerExtModule, only: mem_set_value
@@ -568,7 +625,6 @@ subroutine allocate_scalars(this, name_model, input_mempath)
     character(len=*), intent(in) :: input_mempath
     logical(LGP) :: found
     ! -- local
-! ------------------------------------------------------------------------------
     !
     ! -- Create memory path
     this%memoryPath = create_mem_path(name_model, 'DIS')
@@ -615,21 +671,15 @@ subroutine allocate_scalars(this, name_model, input_mempath)
     !
     ! -- Return
     return
-  end subroutine allocate_scalars
-
-  subroutine allocate_arrays(this)
-! ******************************************************************************
-! allocate_arrays -- Read discretization information from file
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
+  end subroutine allocate_scalars_default
+
+  !> @brief Allocate and initialize array variables.
+  subroutine allocate_arrays_default(this)
     ! -- modules
     use MemoryManagerModule, only: mem_allocate
     ! -- dummy
     class(DisBaseType) :: this
     integer :: isize
-! ------------------------------------------------------------------------------
     !
     ! -- Allocate
     call mem_allocate(this%mshape, this%ndim, 'MSHAPE', this%memoryPath)
@@ -655,88 +705,16 @@ subroutine allocate_arrays(this)
     !
     ! -- Return
     return
-  end subroutine allocate_arrays
-
-  function nodeu_from_string(this, lloc, istart, istop, in, iout, line, &
-                             flag_string, allow_zero) result(nodeu)
-! ******************************************************************************
-! nodeu_from_string -- Receive a string and convert the string to a user
-!   nodenumber.  The model is unstructured; just read user nodenumber.
-!   If flag_string argument is present and true, the first token in string
-!   is allowed to be a string (e.g. boundary name). In this case, if a string
-!   is encountered, return value as -2.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- dummy
-    class(DisBaseType) :: 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
-! ------------------------------------------------------------------------------
-    !
-    !
-    nodeu = 0
-    call store_error('Program error: DisBaseType method nodeu_from_string &
-                     &not implemented.', terminate=.TRUE.)
-    !
-    ! -- return
-    return
-  end function nodeu_from_string
-
-  function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, &
-                             allow_zero) result(nodeu)
-! ******************************************************************************
-! nodeu_from_cellid -- Receive cellid as a string and convert the string to a
-!   user nodenumber.
-!   If flag_string argument is present and true, the first token in string
-!   is allowed to be a string (e.g. boundary name). In this case, if a string
-!   is encountered, return value as -2.
-!   If allow_zero argument is present and true, if all indices equal zero, the
-!   result can be zero. If allow_zero is false, a zero in any index causes an
-!   error.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- dummy
-    class(DisBaseType) :: 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
-    integer(I4B) :: nodeu
-! ------------------------------------------------------------------------------
-    !
-    nodeu = 0
-    call store_error('Program error: DisBaseType method nodeu_from_cellid &
-                      &not implemented.', terminate=.TRUE.)
-    !
-    ! -- return
-    return
-  end function nodeu_from_cellid
+  end subroutine allocate_arrays_default
 
+  !> @brief Parse a reduced nodenumber from a line of text.
+  !!
+  !! The model is unstructured; just read user nodenumber.
+  !! If flag_string argument is present and true, the first token in string
+  !! is allowed to be a string (e.g. boundary name). In this case, if a string
+  !! is encountered, return value as -2.
   function noder_from_string(this, lloc, istart, istop, in, iout, line, &
                              flag_string) result(noder)
-! ******************************************************************************
-! noder_from_string -- Receive a string and convert the string to a reduced
-!   nodenumber.  The model is unstructured; just read user nodenumber.
-!   If flag_string argument is present and true, the first token in string
-!   is allowed to be a string (e.g. boundary name). In this case, if a string
-!   is encountered, return value as -2.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- dummy
     class(DisBaseType) :: this
     integer(I4B), intent(inout) :: lloc
@@ -751,7 +729,6 @@ function noder_from_string(this, lloc, istart, istop, in, iout, line, &
     integer(I4B) :: nodeu
     character(len=LINELENGTH) :: nodestr
     logical :: flag_string_local
-! ------------------------------------------------------------------------------
     !
     if (present(flag_string)) then
       flag_string_local = flag_string
@@ -779,21 +756,15 @@ function noder_from_string(this, lloc, istart, istop, in, iout, line, &
     return
   end function noder_from_string
 
+  !> @brief Parse a reduced nodenumber from a cell ID string.
+  !!
+  !! If flag_string argument is present and true, the first token in string
+  !! is allowed to be a string (e.g. boundary name). In this case, if a string
+  !! is encountered, return value as -2.
+  !! If allow_zero argument is present and true, if all indices equal zero, the
+  !! result can be zero. If allow_zero is false, a zero in any index is an error.
   function noder_from_cellid(this, cellid, inunit, iout, flag_string, &
                              allow_zero) result(noder)
-! ******************************************************************************
-! noder_from_cellid -- Receive cellid as a string and convert it to a reduced
-!   nodenumber.
-!   If flag_string argument is present and true, the first token in string
-!   is allowed to be a string (e.g. boundary name). In this case, if a string
-!   is encountered, return value as -2.
-!   If allow_zero argument is present and true, if all indices equal zero, the
-!   result can be zero. If allow_zero is false, a zero in any index causes an
-!   error.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- return
     integer(I4B) :: noder
     ! -- dummy
@@ -808,7 +779,6 @@ function noder_from_cellid(this, cellid, inunit, iout, flag_string, &
     logical :: allowzerolocal
     character(len=LINELENGTH) :: nodestr
     logical :: flag_string_local
-! ------------------------------------------------------------------------------
     !
     if (present(flag_string)) then
       flag_string_local = flag_string
@@ -842,55 +812,8 @@ function noder_from_cellid(this, cellid, inunit, iout, flag_string, &
     return
   end function noder_from_cellid
 
-  logical function supports_layers(this)
-! ******************************************************************************
-! supports_layers
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- dummy
-    class(DisBaseType) :: this
-! ------------------------------------------------------------------------------
-    !
-    !
-    supports_layers = .false.
-    call store_error('Program error: DisBaseType method supports_layers not &
-                     &implemented.', terminate=.TRUE.)
-    return
-  end function supports_layers
-
-  function get_ncpl(this)
-! ******************************************************************************
-! get_ncpl -- Return number of cells per layer.  This is nodes
-!   for a DISU grid, as there are no layers.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- modules
-    ! -- return
-    integer(I4B) :: get_ncpl
-    ! -- dummy
-    class(DisBaseType) :: this
-! ------------------------------------------------------------------------------
-    !
-    !
-    get_ncpl = 0
-    call store_error('Program error: DisBaseType method get_ncpl not &
-                     &implemented.', terminate=.TRUE.)
-    !
-    ! -- Return
-    return
-  end function get_ncpl
-
+  !> @brief Get the volume of cell n based on x value passed.
   function get_cell_volume(this, n, x)
-! ******************************************************************************
-! get_cell_volume -- Return volume of cell n based on x value passed.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- modules
     ! -- return
     real(DP) :: get_cell_volume
@@ -903,7 +826,6 @@ function get_cell_volume(this, n, x)
     real(DP) :: bt
     real(DP) :: sat
     real(DP) :: thick
-! ------------------------------------------------------------------------------
     !
     get_cell_volume = DZERO
     tp = this%top(n)
@@ -916,14 +838,9 @@ function get_cell_volume(this, n, x)
     return
   end function get_cell_volume
 
+  !> @brief Read an integer array from a line of text.
   subroutine read_int_array(this, line, lloc, istart, istop, iout, in, &
                             iarray, aname)
-! ******************************************************************************
-! read_int_array -- Read a GWF integer array
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- dummy
     class(DisBaseType), intent(inout) :: this
     character(len=*), intent(inout) :: line
@@ -944,14 +861,9 @@ subroutine read_int_array(this, line, lloc, istart, istop, iout, in, &
     return
   end subroutine read_int_array
 
+  !> @brief Read a double precision array from a line of text.
   subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, &
                             darray, aname)
-! ******************************************************************************
-! read_dbl_array -- Read a GWF double precision array
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- dummy
     class(DisBaseType), intent(inout) :: this
     character(len=*), intent(inout) :: line
@@ -972,13 +884,8 @@ subroutine read_dbl_array(this, line, lloc, istart, istop, iout, in, &
     return
   end subroutine read_dbl_array
 
+  !> @brief Fill an integer array indexed by reduced node number.
   subroutine fill_int_array(this, ibuff1, ibuff2)
-! ******************************************************************************
-! fill_dbl_array -- Fill a GWF integer array
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- dummy
     class(DisBaseType), intent(inout) :: this
     integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibuff1
@@ -986,7 +893,7 @@ subroutine fill_int_array(this, ibuff1, ibuff2)
     ! -- local
     integer(I4B) :: nodeu
     integer(I4B) :: noder
-! ------------------------------------------------------------------------------
+    !
     do nodeu = 1, this%nodesuser
       noder = this%get_nodenumber(nodeu, 0)
       if (noder <= 0) cycle
@@ -997,13 +904,8 @@ subroutine fill_int_array(this, ibuff1, ibuff2)
     return
   end subroutine fill_int_array
 
+  !> @brief Fill a double precision array indexed by reduced node number.
   subroutine fill_dbl_array(this, buff1, buff2)
-! ******************************************************************************
-! fill_dbl_array -- Fill a GWF double precision array
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- dummy
     class(DisBaseType), intent(inout) :: this
     real(DP), dimension(:), pointer, contiguous, intent(in) :: buff1
@@ -1011,7 +913,7 @@ subroutine fill_dbl_array(this, buff1, buff2)
     ! -- local
     integer(I4B) :: nodeu
     integer(I4B) :: noder
-! ------------------------------------------------------------------------------
+    !
     do nodeu = 1, this%nodesuser
       noder = this%get_nodenumber(nodeu, 0)
       if (noder <= 0) cycle
@@ -1022,20 +924,16 @@ subroutine fill_dbl_array(this, buff1, buff2)
     return
   end subroutine fill_dbl_array
 
+  !> @brief Read a list using the list reader object.
+  !!
+  !! Convert user node numbers to reduced numbers.
+  !! Terminate if any nodenumbers are within an inactive domain.
+  !! Set up time series and multiply by iauxmultcol if it exists.
+  !! Write the list to iout if iprpak is set.
   subroutine read_list(this, line_reader, in, iout, iprpak, nlist, &
                        inamedbound, iauxmultcol, nodelist, rlist, auxvar, &
                        auxname, boundname, label, pkgname, tsManager, iscloc, &
                        indxconvertflux)
-! ******************************************************************************
-! read_list -- Read a list using the list reader object.
-!              Convert user node numbers to reduced numbers.
-!              Terminate if any nodenumbers are within an inactive domain.
-!              Set up time series and multiply by iauxmultcol if it exists.
-!              Write the list to iout if iprpak is set.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- modules
     use ConstantsModule, only: LENBOUNDNAME, LINELENGTH
     use LongLineReaderModule, only: LongLineReaderType
@@ -1075,7 +973,6 @@ subroutine read_list(this, line_reader, in, iout, iprpak, nlist, &
     type(ListReaderType) :: lstrdobj
     type(TimeSeriesLinkType), pointer :: tsLinkBnd => null()
     type(TimeSeriesLinkType), pointer :: tsLinkAux => null()
-! ------------------------------------------------------------------------------
     !
     ! -- Read the list
     call lstrdobj%read_list(line_reader, in, iout, nlist, inamedbound, &
@@ -1180,86 +1077,8 @@ subroutine read_list(this, line_reader, in, iout, iprpak, nlist, &
     ! -- return
   end subroutine read_list
 
-  subroutine read_layer_array(this, nodelist, darray, ncolbnd, maxbnd, &
-                              icolbnd, aname, inunit, iout)
-! ******************************************************************************
-! read_layer_array -- Read a 2d double array into col icolbnd of darray.
-!                     For cells that are outside of the active domain,
-!                     do not copy the array value into darray.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- modules
-    ! -- dummy
-    class(DisBaseType) :: this
-    integer(I4B), intent(in) :: ncolbnd
-    integer(I4B), intent(in) :: maxbnd
-    integer(I4B), dimension(maxbnd) :: nodelist
-    real(DP), dimension(ncolbnd, maxbnd), intent(inout) :: darray
-    integer(I4B), intent(in) :: icolbnd
-    character(len=*), intent(in) :: aname
-    integer(I4B), intent(in) :: inunit
-    integer(I4B), intent(in) :: iout
-    !
-    !
-    errmsg = 'Programmer error: read_layer_array needs to be overridden &
-            &in any DIS type that extends DisBaseType'
-    call store_error(errmsg, terminate=.TRUE.)
-    !
-    ! -- return
-  end subroutine read_layer_array
-
-  subroutine record_array(this, darray, iout, iprint, idataun, aname, &
-                          cdatafmp, nvaluesp, nwidthp, editdesc, dinact)
-! ******************************************************************************
-! record_array -- Record a double precision array.  The array will be
-!   printed to an external file and/or written to an unformatted external file
-!   depending on the argument specifications.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-!      darray is the double precision array to record
-!      iout is the unit number for ascii output
-!      iprint is a flag indicating whether or not to print the array
-!      idataun is the unit number to which the array will be written in binary
-!        form; if negative then do not write by layers, write entire array
-!      aname is the text descriptor of the array
-!      cdatafmp is the fortran format for writing the array
-!      nvaluesp is the number of values per line for printing
-!      nwidthp is the width of the number for printing
-!      editdesc is the format type (I, G, F, S, E)
-!      dinact is the double precision value to use for cells that are excluded
-!        from the model domain
-! ------------------------------------------------------------------------------
-    ! -- dummy
-    class(DisBaseType), intent(inout) :: this
-    real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray
-    integer(I4B), intent(in) :: iout
-    integer(I4B), intent(in) :: iprint
-    integer(I4B), intent(in) :: idataun
-    character(len=*), intent(in) :: aname
-    character(len=*), intent(in) :: cdatafmp
-    integer(I4B), intent(in) :: nvaluesp
-    integer(I4B), intent(in) :: nwidthp
-    character(len=*), intent(in) :: editdesc
-    real(DP), intent(in) :: dinact
-    !
-    ! --
-    errmsg = 'Programmer error: record_array needs to be overridden &
-            &in any DIS type that extends DisBaseType'
-    call store_error(errmsg, terminate=.TRUE.)
-    !
-  end subroutine record_array
-
+  !> @brief Record a connection-based double precision array.
   subroutine record_connection_array(this, flowja, ibinun, iout)
-! ******************************************************************************
-! record_connection_array -- Record a connection-based double precision
-! array for either a structured or an unstructured grid.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- dummy
     class(DisBaseType) :: this
     real(DP), dimension(:), intent(in) :: flowja
@@ -1269,7 +1088,6 @@ subroutine record_connection_array(this, flowja, ibinun, iout)
     character(len=16), dimension(1) :: text
     ! -- data
     data text(1)/'    FLOW-JA-FACE'/
-! ------------------------------------------------------------------------------
     !
     ! -- write full ja array
     call ubdsv1(kstp, kper, text(1), ibinun, flowja, size(flowja), 1, 1, &
@@ -1279,14 +1097,8 @@ subroutine record_connection_array(this, flowja, ibinun, iout)
     return
   end subroutine record_connection_array
 
+  !> @brief Convert reduced node number to string in form (nodenumber) or (k,i,j).
   subroutine noder_to_string(this, noder, str)
-! ******************************************************************************
-! noder_to_string -- Convert reduced node number to a string in the form of
-! (nodenumber) or (k,i,j)
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- modules
     ! -- dummy
     class(DisBaseType) :: this
@@ -1294,7 +1106,6 @@ subroutine noder_to_string(this, noder, str)
     character(len=*), intent(inout) :: str
     ! -- local
     integer(I4B) :: nodeu
-! ------------------------------------------------------------------------------
     !
     nodeu = this%get_nodeuser(noder)
     call this%nodeu_to_string(nodeu, str)
@@ -1303,14 +1114,9 @@ subroutine noder_to_string(this, noder, str)
     return
   end subroutine noder_to_string
 
+  !> @brief Convert reduced node number to cellid and fill array with (nodenumber)
+  !! or (k,j) or (k,i,j)
   subroutine noder_to_array(this, noder, arr)
-! ******************************************************************************
-! noder_to_array -- Convert reduced node number to cellid and fill array with
-!                   (nodenumber) or (k,j) or (k,i,j)
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- modules
     ! -- dummy
     class(DisBaseType) :: this
@@ -1318,7 +1124,6 @@ subroutine noder_to_array(this, noder, arr)
     integer(I4B), dimension(:), intent(inout) :: arr
     ! -- local
     integer(I4B) :: nodeu
-! ------------------------------------------------------------------------------
     !
     nodeu = this%get_nodeuser(noder)
     call this%nodeu_to_array(nodeu, arr)
@@ -1327,45 +1132,9 @@ subroutine noder_to_array(this, noder, arr)
     return
   end subroutine noder_to_array
 
-  subroutine record_srcdst_list_header(this, text, textmodel, textpackage, &
-                                       dstmodel, dstpackage, naux, auxtxt, &
-                                       ibdchn, nlist, iout)
-! ******************************************************************************
-! record_srcdst_list_header -- Record list header for imeth=6
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- dummy
-    class(DisBaseType) :: 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
-    !
-    ! --
-    errmsg = 'Programmer error: record_srcdst_list_header needs to be &
-            &overridden in any DIS type that extends DisBaseType'
-    call store_error(errmsg, terminate=.TRUE.)
-    !
-    ! -- return
-    return
-  end subroutine record_srcdst_list_header
-
+  !> @brief Record list header
   subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, &
                                       naux, aux, olconv, olconv2)
-! ******************************************************************************
-! record_srcdst_list_header -- Record list header
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- modules
     use InputOutputModule, only: ubdsvd
     ! -- dummy
@@ -1383,7 +1152,6 @@ subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, &
     logical :: lconv2
     integer(I4B) :: nodeu
     integer(I4B) :: nodeu2
-! ------------------------------------------------------------------------------
     !
     ! -- Use ubdsvb to write list header
     if (present(olconv)) then
@@ -1412,42 +1180,8 @@ subroutine record_srcdst_list_entry(this, ibdchn, noder, noder2, q, &
     return
   end subroutine record_srcdst_list_entry
 
-  subroutine nlarray_to_nodelist(this, darray, nodelist, maxbnd, nbound, aname)
-! ******************************************************************************
-! nlarray_to_nodelist -- Convert an integer array into nodelist. For structured
-!                        model, integer array is layer number; for unstructured
-!                        model, integer array is node number.
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    ! -- modules
-    use SimModule, only: store_error
-    use ConstantsModule, only: LINELENGTH
-    ! -- dummy
-    class(DisBaseType) :: this
-    integer(I4B), intent(in) :: maxbnd
-    integer(I4B), dimension(:), pointer, contiguous :: darray
-    integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
-    integer(I4B), intent(inout) :: nbound
-    character(len=*), intent(in) :: aname
-    !
-    ! --
-    errmsg = 'Programmer error: nlarray_to_nodelist needs to be &
-            &overridden in any DIS type that extends DisBaseType'
-    call store_error(errmsg, terminate=.TRUE.)
-    !
-    ! -- return
-    return
-  end subroutine nlarray_to_nodelist
-
+  !> @brief Find the first highest active cell beneath cell n.
   subroutine highest_active(this, n, ibound)
-! ******************************************************************************
-! highest_active -- Find the first highest active cell beneath cell n
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- dummy
     class(DisBaseType) :: this
     integer(I4B), intent(inout) :: n
@@ -1455,7 +1189,6 @@ subroutine highest_active(this, n, ibound)
     ! -- locals
     integer(I4B) :: m, ii, iis
     logical done, bottomcell
-! ------------------------------------------------------------------------------
     !
     ! -- Loop through connected cells until the highest active one (including a
     !    constant head cell) is found.  Return that cell as n.
@@ -1488,19 +1221,13 @@ subroutine highest_active(this, n, ibound)
     return
   end subroutine highest_active
 
+  !> @brief Get the cell area for this node.
   function get_area(this, node) result(area)
-! ******************************************************************************
-! get_area -- Return the cell area for this node
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
     ! -- return
     real(DP) :: area
     ! -- dummy
     class(DisBaseType) :: this
     integer(I4B), intent(in) :: node
-! ------------------------------------------------------------------------------
     !
     ! -- Return the cell area
     area = this%area(node)
@@ -1509,16 +1236,14 @@ function get_area(this, node) result(area)
     return
   end function get_area
 
-  !> @ brief Calculate the area factor for the cell connection
-  !!
-  !!  Function calculates the area factor for the cell connection. The sum of
-  !!  all area factors for all cell connections to overlying or underlying
-  !!  cells cells will be 1.
+  !> @brief Calculate the area factor for the cell connection
   !!
-  !!  TODO: confirm that this works for cells that are only partially covered
-  !!        by overlying or underlying cells.
+  !! Function calculates the area factor for the cell connection. The sum of
+  !! all area factors for all cell connections to overlying or underlying
+  !! cells cells will be 1.
   !!
-  !<
+  !! TODO: confirm that this works for cells that are only partially covered
+  !!       by overlying or underlying cells.
   function get_area_factor(this, node, idx_conn) result(area_factor)
     ! -- return
     real(DP) :: area_factor !< connection cell area factor
@@ -1529,7 +1254,6 @@ function get_area_factor(this, node, idx_conn) result(area_factor)
     ! -- local
     real(DP) :: area_node
     real(DP) :: area_conn
-    ! ------------------------------------------------------------------------------
     !
     ! -- calculate the cell area fraction
     area_node = this%area(node)