diff --git a/src/Exchange/exg-gwfgwe.f90 b/src/Exchange/exg-gwfgwe.f90
index 536a57b232f..485a6bc5b82 100644
--- a/src/Exchange/exg-gwfgwe.f90
+++ b/src/Exchange/exg-gwfgwe.f90
@@ -229,32 +229,10 @@ subroutine exg_ar(this)
     end if
     !
     ! -- Make sure idomains are identical
-    select type (gwfdis => gwfmodel%dis)
-    type is (DisType)
-      select type (gwedis => gwemodel%dis)
-      type is (DisType)
-        if (.not. all(gwfdis%idomain == gwedis%idomain)) then
-          write (errmsg, fmtidomerr) trim(this%name)
-          call store_error(errmsg, terminate=.TRUE.)
-        end if
-      end select
-    type is (DisvType)
-      select type (gwedis => gwemodel%dis)
-      type is (DisvType)
-        if (.not. all(gwfdis%idomain == gwedis%idomain)) then
-          write (errmsg, fmtidomerr) trim(this%name)
-          call store_error(errmsg, terminate=.TRUE.)
-        end if
-      end select
-    type is (DisuType)
-      select type (gwedis => gwemodel%dis)
-      type is (DisuType)
-        if (.not. all(gwfdis%idomain == gwedis%idomain)) then
-          write (errmsg, fmtidomerr) trim(this%name)
-          call store_error(errmsg, terminate=.TRUE.)
-        end if
-      end select
-    end select
+    if (.not. all(gwfmodel%dis%idomain == gwemodel%dis%idomain)) then
+      write (errmsg, fmtidomerr) trim(this%name)
+      call store_error(errmsg, terminate=.TRUE.)
+    end if
     !
     ! -- setup pointers to gwf variables allocated in gwf_ar
     gwemodel%fmi%gwfhead => gwfmodel%x
diff --git a/src/Exchange/exg-gwfgwt.f90 b/src/Exchange/exg-gwfgwt.f90
index aa198e7a100..67f29a26672 100644
--- a/src/Exchange/exg-gwfgwt.f90
+++ b/src/Exchange/exg-gwfgwt.f90
@@ -232,32 +232,10 @@ subroutine exg_ar(this)
     end if
     !
     ! -- Make sure idomains are identical
-    select type (gwfdis => gwfmodel%dis)
-    type is (DisType)
-      select type (gwtdis => gwtmodel%dis)
-      type is (DisType)
-        if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
-          write (errmsg, fmtidomerr) trim(this%name)
-          call store_error(errmsg, terminate=.TRUE.)
-        end if
-      end select
-    type is (DisvType)
-      select type (gwtdis => gwtmodel%dis)
-      type is (DisvType)
-        if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
-          write (errmsg, fmtidomerr) trim(this%name)
-          call store_error(errmsg, terminate=.TRUE.)
-        end if
-      end select
-    type is (DisuType)
-      select type (gwtdis => gwtmodel%dis)
-      type is (DisuType)
-        if (.not. all(gwfdis%idomain == gwtdis%idomain)) then
-          write (errmsg, fmtidomerr) trim(this%name)
-          call store_error(errmsg, terminate=.TRUE.)
-        end if
-      end select
-    end select
+    if (.not. all(gwfmodel%dis%idomain == gwtmodel%dis%idomain)) then
+      write (errmsg, fmtidomerr) trim(this%name)
+      call store_error(errmsg, terminate=.TRUE.)
+    end if
     !
     ! -- setup pointers to gwf variables allocated in gwf_ar
     gwtmodel%fmi%gwfhead => gwfmodel%x
diff --git a/src/Exchange/exg-gwfprt.f90 b/src/Exchange/exg-gwfprt.f90
index 61929f13693..f60cd58f582 100644
--- a/src/Exchange/exg-gwfprt.f90
+++ b/src/Exchange/exg-gwfprt.f90
@@ -225,32 +225,10 @@ subroutine exg_ar(this)
     end if
     !
     ! -- Make sure idomains are identical
-    select type (gwfdis => gwfmodel%dis)
-    type is (DisType)
-      select type (prtdis => prtmodel%dis)
-      type is (DisType)
-        if (.not. all(gwfdis%idomain == prtdis%idomain)) then
-          write (errmsg, fmtidomerr) trim(this%name)
-          call store_error(errmsg, terminate=.TRUE.)
-        end if
-      end select
-    type is (DisvType)
-      select type (prtdis => prtmodel%dis)
-      type is (DisvType)
-        if (.not. all(gwfdis%idomain == prtdis%idomain)) then
-          write (errmsg, fmtidomerr) trim(this%name)
-          call store_error(errmsg, terminate=.TRUE.)
-        end if
-      end select
-    type is (DisuType)
-      select type (prtdis => prtmodel%dis)
-      type is (DisuType)
-        if (.not. all(gwfdis%idomain == prtdis%idomain)) then
-          write (errmsg, fmtidomerr) trim(this%name)
-          call store_error(errmsg, terminate=.TRUE.)
-        end if
-      end select
-    end select
+    if (.not. all(gwfmodel%dis%idomain == prtmodel%dis%idomain)) then
+      write (errmsg, fmtidomerr) trim(this%name)
+      call store_error(errmsg, terminate=.TRUE.)
+    end if
     !
     ! -- setup pointers to gwf variables allocated in gwf_ar
     prtmodel%fmi%gwfhead => gwfmodel%x
diff --git a/src/Model/Discretization/Dis.f90 b/src/Model/Discretization/Dis.f90
index 798d1ae479f..a88b99b02b9 100644
--- a/src/Model/Discretization/Dis.f90
+++ b/src/Model/Discretization/Dis.f90
@@ -28,11 +28,9 @@ module DisModule
     real(DP), dimension(:), pointer, contiguous :: delc => null() !< spacing along a column
     real(DP), dimension(:, :), pointer, contiguous :: top2d => null() !< top elevations for each cell at top of model (ncol, nrow)
     real(DP), dimension(:, :, :), pointer, contiguous :: bot3d => null() !< bottom elevations for each cell (ncol, nrow, nlay)
-    integer(I4B), dimension(:, :, :), pointer, contiguous :: idomain => null() !< idomain (ncol, nrow, nlay)
     real(DP), dimension(:, :, :), pointer :: botm => null() !< top and bottom elevations for each cell (ncol, nrow, nlay)
     real(DP), dimension(:), pointer, contiguous :: cellx => null() !< cell center x coordinate for column j
     real(DP), dimension(:), pointer, contiguous :: celly => null() !< cell center y coordinate for row i
-
   contains
 
     procedure :: dis_df => dis3d_df
@@ -244,7 +242,6 @@ subroutine source_dimensions(this)
     ! -- dummy
     class(DisType) :: this
     ! -- locals
-    integer(I4B) :: i, j, k
     type(DisFoundType) :: found
     !
     ! -- update defaults with idm sourced values
@@ -280,23 +277,13 @@ subroutine source_dimensions(this)
     ! -- Allocate delr, delc, and non-reduced vectors for dis
     call mem_allocate(this%delr, this%ncol, 'DELR', this%memoryPath)
     call mem_allocate(this%delc, this%nrow, 'DELC', this%memoryPath)
-    call mem_allocate(this%idomain, this%ncol, this%nrow, this%nlay, 'IDOMAIN', &
-                      this%memoryPath)
+    call mem_allocate(this%idomain, this%nodesuser, 'IDOMAIN', this%memoryPath)
     call mem_allocate(this%top2d, this%ncol, this%nrow, 'TOP2D', this%memoryPath)
     call mem_allocate(this%bot3d, this%ncol, this%nrow, this%nlay, 'BOT3D', &
                       this%memoryPath)
     call mem_allocate(this%cellx, this%ncol, 'CELLX', this%memoryPath)
     call mem_allocate(this%celly, this%nrow, 'CELLY', this%memoryPath)
     !
-    ! -- initialize all cells to be active (idomain = 1)
-    do k = 1, this%nlay
-      do i = 1, this%nrow
-        do j = 1, this%ncol
-          this%idomain(j, i, k) = 1
-        end do
-      end do
-    end do
-    !
   end subroutine source_dimensions
 
   !> @brief Write dimensions to list file
@@ -329,14 +316,28 @@ end subroutine log_dimensions
   subroutine source_griddata(this)
     ! -- dummy
     class(DisType) :: this
+    ! -- local
     type(DisFoundType) :: found
+    integer(I4B) :: j, i, k
+    integer(I4B), contiguous, pointer :: idomain(:, :, :)
+    !
+    ! -- initialize all cells to be active (idomain = 1)
+    allocate (idomain(this%ncol, this%nrow, this%nlay))
+    do k = 1, this%nlay
+      do i = 1, this%nrow
+        do j = 1, this%ncol
+          idomain(j, i, k) = 1
+        end do
+      end do
+    end do
     !
     ! -- update defaults with idm sourced values
     call mem_set_value(this%delr, 'DELR', this%input_mempath, found%delr)
     call mem_set_value(this%delc, 'DELC', this%input_mempath, found%delc)
     call mem_set_value(this%top2d, 'TOP', this%input_mempath, found%top)
     call mem_set_value(this%bot3d, 'BOTM', this%input_mempath, found%botm)
-    call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    call mem_set_value(idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    this%idomain = reshape(idomain, [this%nodesuser])
     !
     ! -- log simulation values
     if (this%iout > 0) then
@@ -386,7 +387,7 @@ subroutine grid_finalize(this)
     ! -- dummy
     class(DisType) :: this
     ! -- locals
-    integer(I4B) :: n, i, j, k
+    integer(I4B) :: n, nn, k, i, j
     integer(I4B) :: node
     integer(I4B) :: noder
     integer(I4B) :: nrsize
@@ -403,12 +404,8 @@ subroutine grid_finalize(this)
     !
     ! -- count active cells
     this%nodes = 0
-    do k = 1, this%nlay
-      do i = 1, this%nrow
-        do j = 1, this%ncol
-          if (this%idomain(j, i, k) > 0) this%nodes = this%nodes + 1
-        end do
-      end do
+    do i = 1, this%nodesuser
+      if (this%idomain(i) > 0) this%nodes = this%nodes + 1
     end do
     !
     ! -- Check to make sure nodes is a valid number
@@ -424,7 +421,8 @@ subroutine grid_finalize(this)
     do k = 1, this%nlay
       do i = 1, this%nrow
         do j = 1, this%ncol
-          if (this%idomain(j, i, k) < 1) cycle
+          nn = get_node(k, i, j, this%nlay, this%nrow, this%ncol)
+          if (this%idomain(nn) < 1) cycle
           if (k > 1) then
             top = this%bot3d(j, i, k - 1)
           else
@@ -461,10 +459,11 @@ subroutine grid_finalize(this)
       do k = 1, this%nlay
         do i = 1, this%nrow
           do j = 1, this%ncol
-            if (this%idomain(j, i, k) > 0) then
+            nn = get_node(k, i, j, this%nlay, this%nrow, this%ncol)
+            if (this%idomain(nn) > 0) then
               this%nodereduced(node) = noder
               noder = noder + 1
-            elseif (this%idomain(j, i, k) < 0) then
+            elseif (this%idomain(nn) < 0) then
               this%nodereduced(node) = -1
             else
               this%nodereduced(node) = 0
@@ -482,7 +481,8 @@ subroutine grid_finalize(this)
       do k = 1, this%nlay
         do i = 1, this%nrow
           do j = 1, this%ncol
-            if (this%idomain(j, i, k) > 0) then
+            nn = get_node(k, i, j, this%nlay, this%nrow, this%ncol)
+            if (this%idomain(nn) > 0) then
               this%nodeuser(noder) = node
               noder = noder + 1
             end if
diff --git a/src/Model/Discretization/Dis2d.f90 b/src/Model/Discretization/Dis2d.f90
index 6e1446e2c8a..4d7c8b3517f 100644
--- a/src/Model/Discretization/Dis2d.f90
+++ b/src/Model/Discretization/Dis2d.f90
@@ -26,7 +26,6 @@ module Dis2dModule
     real(DP), dimension(:), pointer, contiguous :: delr => null() !< spacing along a row
     real(DP), dimension(:), pointer, contiguous :: delc => null() !< spacing along a column
     real(DP), dimension(:, :), pointer, contiguous :: bottom => null() !< bottom elevations for each cell (ncol, nrow)
-    integer(I4B), dimension(:, :), pointer, contiguous :: idomain => null() !< idomain (ncol, nrow)
     real(DP), dimension(:), pointer, contiguous :: cellx => null() !< cell center x coordinate for column j
     real(DP), dimension(:), pointer, contiguous :: celly => null() !< cell center y coordinate for row i
 
@@ -237,7 +236,6 @@ subroutine source_dimensions(this)
     ! -- dummy
     class(Dis2dType) :: this
     ! -- locals
-    integer(I4B) :: i, j
     type(DisFoundType) :: found
     !
     ! -- update defaults with idm sourced values
@@ -267,20 +265,12 @@ subroutine source_dimensions(this)
     ! -- Allocate delr, delc, and non-reduced vectors for dis
     call mem_allocate(this%delr, this%ncol, 'DELR', this%memoryPath)
     call mem_allocate(this%delc, this%nrow, 'DELC', this%memoryPath)
-    call mem_allocate(this%idomain, this%ncol, this%nrow, 'IDOMAIN', &
-                      this%memoryPath)
+    call mem_allocate(this%idomain, this%nodesuser, 'IDOMAIN', this%memoryPath)
     call mem_allocate(this%bottom, this%ncol, this%nrow, 'BOTTOM', &
                       this%memoryPath)
     call mem_allocate(this%cellx, this%ncol, 'CELLX', this%memoryPath)
     call mem_allocate(this%celly, this%nrow, 'CELLY', this%memoryPath)
     !
-    ! -- initialize all cells to be active (idomain = 1)
-    do i = 1, this%nrow
-      do j = 1, this%ncol
-        this%idomain(j, i) = 1
-      end do
-    end do
-    !
   end subroutine source_dimensions
 
   !> @brief Write dimensions to list file
@@ -309,13 +299,25 @@ end subroutine log_dimensions
   subroutine source_griddata(this)
     ! -- dummy
     class(Dis2dType) :: this
+    ! -- local
     type(DisFoundType) :: found
+    integer(I4B) :: i, j
+    integer(I4B), contiguous, pointer :: idomain(:, :)
+    !
+    ! -- initialize all cells to be active (idomain = 1)
+    allocate (idomain(this%ncol, this%nrow))
+    do i = 1, this%nrow
+      do j = 1, this%ncol
+        idomain(j, i) = 1
+      end do
+    end do
     !
     ! -- update defaults with idm sourced values
     call mem_set_value(this%delr, 'DELR', this%input_mempath, found%delr)
     call mem_set_value(this%delc, 'DELC', this%input_mempath, found%delc)
     call mem_set_value(this%bottom, 'BOTTOM', this%input_mempath, found%bottom)
-    call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    call mem_set_value(idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    this%idomain = reshape(idomain, [this%nodesuser])
     !
     ! -- log simulation values
     if (this%iout > 0) then
@@ -376,10 +378,8 @@ subroutine grid_finalize(this)
     !
     ! -- count active cells
     this%nodes = 0
-    do i = 1, this%nrow
-      do j = 1, this%ncol
-        if (this%idomain(j, i) > 0) this%nodes = this%nodes + 1
-      end do
+    do i = 1, this%nodesuser
+      if (this%idomain(i) > 0) this%nodes = this%nodes + 1
     end do
     !
     ! -- Check to make sure nodes is a valid number
@@ -407,10 +407,10 @@ subroutine grid_finalize(this)
       noder = 1
       do i = 1, this%nrow
         do j = 1, this%ncol
-          if (this%idomain(j, i) > 0) then
+          if (this%idomain(node) > 0) then
             this%nodereduced(node) = noder
             noder = noder + 1
-          elseif (this%idomain(j, i) < 0) then
+          elseif (this%idomain(node) < 0) then
             this%nodereduced(node) = -1
           else
             this%nodereduced(node) = 0
@@ -426,7 +426,7 @@ subroutine grid_finalize(this)
       noder = 1
       do i = 1, this%nrow
         do j = 1, this%ncol
-          if (this%idomain(j, i) > 0) then
+          if (this%idomain(node) > 0) then
             this%nodeuser(noder) = node
             noder = noder + 1
           end if
diff --git a/src/Model/Discretization/Disu.f90 b/src/Model/Discretization/Disu.f90
index 15b2df64058..d878e596676 100644
--- a/src/Model/Discretization/Disu.f90
+++ b/src/Model/Discretization/Disu.f90
@@ -43,9 +43,7 @@ module DisuModule
     integer(I4B), pointer :: iangledegx => null() ! =1 when angle information was present in input, 0 otherwise
     integer(I4B), dimension(:), pointer, contiguous :: iavert => null() ! cell vertex pointer ia array
     integer(I4B), dimension(:), pointer, contiguous :: javert => null() ! cell vertex pointer ja array
-    integer(I4B), dimension(:), pointer, contiguous :: idomain => null() ! idomain (nodes)
     logical(LGP) :: readFromFile ! True, when DIS is read from file (almost always)
-
   contains
 
     procedure :: dis_df => disu_df
diff --git a/src/Model/Discretization/Disv.f90 b/src/Model/Discretization/Disv.f90
index 6781a156175..01643484e09 100644
--- a/src/Model/Discretization/Disv.f90
+++ b/src/Model/Discretization/Disv.f90
@@ -31,8 +31,6 @@ module DisvModule
     integer(I4B), dimension(:), pointer, contiguous :: javert => null() !< cell vertex pointer ja array
     real(DP), dimension(:), pointer, contiguous :: top1d => null() !< top elevations for each cell at top of model (ncpl)
     real(DP), dimension(:, :), pointer, contiguous :: bot2d => null() !< bottom elevations for each cell (ncpl, nlay)
-    integer(I4B), dimension(:, :), pointer, contiguous :: idomain => null() !< idomain (ncpl, nlay)
-
   contains
 
     procedure :: dis_df => disv_df
@@ -262,7 +260,7 @@ subroutine source_dimensions(this)
     ! -- dummy
     class(DisvType) :: this
     ! -- locals
-    integer(I4B) :: j, k
+    integer(I4B) :: i
     type(DisvFoundType) :: found
     !
     ! -- update defaults with idm sourced values
@@ -296,8 +294,7 @@ subroutine source_dimensions(this)
     this%nodesuser = this%nlay * this%ncpl
     !
     ! -- Allocate non-reduced vectors for disv
-    call mem_allocate(this%idomain, this%ncpl, this%nlay, 'IDOMAIN', &
-                      this%memoryPath)
+    call mem_allocate(this%idomain, this%nodesuser, 'IDOMAIN', this%memoryPath)
     call mem_allocate(this%top1d, this%ncpl, 'TOP1D', this%memoryPath)
     call mem_allocate(this%bot2d, this%ncpl, this%nlay, 'BOT2D', &
                       this%memoryPath)
@@ -307,10 +304,8 @@ subroutine source_dimensions(this)
     call mem_allocate(this%cellxy, 2, this%ncpl, 'CELLXY', this%memoryPath)
     !
     ! -- initialize all cells to be active (idomain = 1)
-    do k = 1, this%nlay
-      do j = 1, this%ncpl
-        this%idomain(j, k) = 1
-      end do
+    do i = 1, this%nodesuser
+      this%idomain(i) = 1
     end do
     !
   end subroutine source_dimensions
@@ -347,11 +342,22 @@ subroutine source_griddata(this)
     class(DisvType) :: this
     ! -- locals
     type(DisvFoundType) :: found
+    integer(I4B) :: j, k
+    integer(I4B), contiguous, pointer :: idomain(:, :)
+    !
+    ! -- initialize all cells to be active (idomain = 1)
+    allocate (idomain(this%ncpl, this%nlay))
+    do k = 1, this%nlay
+      do j = 1, this%ncpl
+        idomain(j, k) = 1
+      end do
+    end do
     !
     ! -- update defaults with idm sourced values
     call mem_set_value(this%top1d, 'TOP', this%input_mempath, found%top)
     call mem_set_value(this%bot2d, 'BOTM', this%input_mempath, found%botm)
-    call mem_set_value(this%idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    call mem_set_value(idomain, 'IDOMAIN', this%input_mempath, found%idomain)
+    this%idomain = reshape(idomain, [this%nodesuser])
     !
     ! -- log simulation values
     if (this%iout > 0) then
@@ -405,10 +411,8 @@ subroutine grid_finalize(this)
     !
     ! -- count active cells
     this%nodes = 0
-    do k = 1, this%nlay
-      do j = 1, this%ncpl
-        if (this%idomain(j, k) > 0) this%nodes = this%nodes + 1
-      end do
+    do node = 1, this%nodesuser
+      if (this%idomain(node) > 0) this%nodes = this%nodes + 1
     end do
     !
     ! -- Check to make sure nodes is a valid number
@@ -420,10 +424,14 @@ subroutine grid_finalize(this)
     end if
     !
     ! -- Check cell thicknesses
+    node = 1
     do k = 1, this%nlay
       do j = 1, this%ncpl
-        if (this%idomain(j, k) == 0) cycle
-        if (this%idomain(j, k) > 0) then
+        if (this%idomain(node) == 0) then
+          node = node + 1
+          cycle
+        end if
+        if (this%idomain(node) > 0) then
           if (k > 1) then
             top = this%bot2d(j, k - 1)
           else
@@ -435,6 +443,7 @@ subroutine grid_finalize(this)
             call store_error(errmsg)
           end if
         end if
+        node = node + 1
       end do
     end do
     if (count_errors() > 0) then
@@ -458,10 +467,10 @@ subroutine grid_finalize(this)
       noder = 1
       do k = 1, this%nlay
         do j = 1, this%ncpl
-          if (this%idomain(j, k) > 0) then
+          if (this%idomain(node) > 0) then
             this%nodereduced(node) = noder
             noder = noder + 1
-          elseif (this%idomain(j, k) < 0) then
+          elseif (this%idomain(node) < 0) then
             this%nodereduced(node) = -1
           else
             this%nodereduced(node) = 0
@@ -477,7 +486,7 @@ subroutine grid_finalize(this)
       noder = 1
       do k = 1, this%nlay
         do j = 1, this%ncpl
-          if (this%idomain(j, k) > 0) then
+          if (this%idomain(node) > 0) then
             this%nodeuser(noder) = node
             noder = noder + 1
           end if
diff --git a/src/Model/Discretization/Disv1d.f90 b/src/Model/Discretization/Disv1d.f90
index 16ec85d6f1a..6edc301b79a 100644
--- a/src/Model/Discretization/Disv1d.f90
+++ b/src/Model/Discretization/Disv1d.f90
@@ -23,7 +23,6 @@ module Disv1dModule
     real(DP), dimension(:), pointer, contiguous :: length => null() !< length of each reach (of size nodesuser)
     real(DP), dimension(:), pointer, contiguous :: width => null() !< reach width
     real(DP), dimension(:), pointer, contiguous :: bottom => null() !< reach bottom elevation
-    integer(I4B), dimension(:), pointer, contiguous :: idomain => null() !< idomain (of size nodesuser)
     real(DP), dimension(:, :), pointer, contiguous :: vertices => null() !< cell vertices stored as 2d array with columns of x, y
     real(DP), dimension(:, :), pointer, contiguous :: cellxy => null() !< reach midpoints stored as 2d array with columns of x, y
     real(DP), dimension(:), pointer, contiguous :: fdc => null() !< fdc stored as array
diff --git a/src/Model/Discretization/Disv2d.f90 b/src/Model/Discretization/Disv2d.f90
index e0eb136bf27..d11ee872e20 100644
--- a/src/Model/Discretization/Disv2d.f90
+++ b/src/Model/Discretization/Disv2d.f90
@@ -28,8 +28,6 @@ module Disv2dModule
     integer(I4B), dimension(:), pointer, contiguous :: iavert => null() !< cell vertex pointer ia array
     integer(I4B), dimension(:), pointer, contiguous :: javert => null() !< cell vertex pointer ja array
     real(DP), dimension(:), pointer, contiguous :: bottom => null() !< bottom elevations for each cell (nodes)
-    integer(I4B), dimension(:), pointer, contiguous :: idomain => null() !< idomain (nodes)
-
   contains
 
     procedure :: dis_df => disv2d_df
diff --git a/src/Model/ModelUtilities/DiscretizationBase.f90 b/src/Model/ModelUtilities/DiscretizationBase.f90
index 05151ad9dbd..67b22af1d0e 100644
--- a/src/Model/ModelUtilities/DiscretizationBase.f90
+++ b/src/Model/ModelUtilities/DiscretizationBase.f90
@@ -53,6 +53,7 @@ module BaseDisModule
     integer(I4B), dimension(:), pointer, contiguous :: ibuff => null() !< helper int array of size nodesuser
     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)
+    integer(I4B), dimension(:), pointer, contiguous :: idomain => null() !< (size:nodesuser)
   contains
     procedure :: dis_df
     procedure :: dis_ac
diff --git a/utils/mf5to6/src/Preproc/Discretization3D.f90 b/utils/mf5to6/src/Preproc/Discretization3D.f90
index 078ef2ea263..52188052231 100644
--- a/utils/mf5to6/src/Preproc/Discretization3D.f90
+++ b/utils/mf5to6/src/Preproc/Discretization3D.f90
@@ -24,7 +24,7 @@ module DnmDis3dModule
     double precision                            :: gridYmin = DZERO
     double precision, dimension(:), allocatable :: delr
     double precision, dimension(:), allocatable :: delc
-    integer,      dimension(:,:,:), allocatable :: idomain
+    integer,      dimension(:), allocatable :: idomain
   contains
     procedure :: dis_df => dis3d_df
     procedure :: dis_da => dis3d_da
@@ -293,7 +293,7 @@ subroutine read_data(this)
 !    double precision, dimension(:),     allocatable  :: delr
 !    double precision, dimension(:),     allocatable  :: delc
     double precision, dimension(:,:,:), allocatable  :: botm
-!    integer,          dimension(:,:,:), allocatable  :: idomain
+    integer,          dimension(:,:,:), allocatable  :: idomain
     character(len=LINELENGTH) :: keyword
     integer :: lloc, istart, istop, n, node, noder, i, j, k
     integer :: ierr, ival
@@ -324,9 +324,10 @@ subroutine read_data(this)
     ! -- Allocate temporary arrays used in this subroutine
     allocate(this%delr(this%ncol))
     allocate(this%delc(this%nrow))
-    allocate(this%idomain(this%ncol, this%nrow, this%nlay))
-    allocate(botm(this%ncol, this%nrow, 0:this%nlay))
     this%nodesuser = this%nlay * this%nrow * this%ncol
+    allocate(idomain(this%ncol, this%nrow, this%nlay))
+    allocate(this%idomain(this%nodesuser))
+    allocate(botm(this%ncol, this%nrow, 0:this%nlay))
     !
     ! --Read GRIDDATA block
     call this%parser%GetBlock('GRIDDATA', isfound, ierr, supportOpenClose=.true.)
@@ -363,10 +364,11 @@ subroutine read_data(this)
           call this%parser%GetStringCaps(keyword)
           if (keyword.EQ.'LAYERED') then
             call ReadArray(this%inunit, this%iout, this%nlay, this%nrow,        &
-                           this%ncol, this%nodesuser, this%idomain(:, :, :),         &
+                           this%ncol, this%nodesuser, idomain(:, :, :),         &
                            aname(5))
+            this%idomain = reshape(idomain, [this%nodesuser])
           else
-            call ReadArray(this%idomain(:, :, :), aname(5), this%nodesuser,          &
+            call ReadArray(idomain, aname(5), this%nodesuser,          &
                         this%inunit, this%iout)
           end if
           lname(5) = .true.
@@ -401,22 +403,14 @@ subroutine read_data(this)
     ! -- If IDOMAIN was not read, then set all values to 1, otherwise
     !    count active cells
     if(.not. lname(5)) then
-      do k = 1, this%nlay
-        do i = 1, this%nrow
-          do j = 1, this%ncol
-            this%idomain(j, i, k) = 1
-          enddo
-        enddo
+      do node = 1, this%nodesuser
+        this%idomain(node) = 1
       enddo
       this%nodes = this%nodesuser
     else
       this%nodes = 0
-      do k = 1, this%nlay
-        do i = 1, this%nrow
-          do j = 1, this%ncol
-            if(this%idomain(j, i, k) > 0) this%nodes = this%nodes + 1
-          enddo
-        enddo
+      do node = 1, this%nodesuser
+        if(this%idomain(node) > 0) this%nodes = this%nodes + 1
       enddo
     endif
     !
@@ -424,8 +418,9 @@ subroutine read_data(this)
     do k = 1, this%nlay
       do i = 1, this%nrow
         do j = 1, this%ncol
-          if (this%idomain(j, i, k) == 0) cycle
-          if (this%idomain(j, i, k) > 0) then
+          node = get_node(k, i, j, this%nlay, this%nrow, this%ncol)
+          if (this%idomain(node) == 0) cycle
+          if (this%idomain(node) > 0) then
             dz = botm(j, i, k - 1) - botm(j, i, k)
             if (dz <= DZERO) then
               write(ermsg, fmt=fmtdz) k, i, j, botm(j, i, k - 1), botm(j, i, k)
@@ -455,10 +450,10 @@ subroutine read_data(this)
       do k = 1, this%nlay
         do i = 1, this%nrow
           do j = 1, this%ncol
-            if(this%idomain(j, i, k) > 0) then
+            if(this%idomain(node) > 0) then
               this%nodereduced(node) = noder
               noder = noder + 1
-            elseif(this%idomain(j, i, k) < 0) then
+            elseif(this%idomain(node) < 0) then
               this%nodereduced(node) = -1
             else
               this%nodereduced(node) = 0
@@ -476,7 +471,7 @@ subroutine read_data(this)
       do k = 1, this%nlay
         do i = 1, this%nrow
           do j = 1, this%ncol
-            if(this%idomain(j, i, k) > 0) then
+            if(this%idomain(node) > 0) then
               this%nodeuser(noder) = node
               noder = noder + 1
             endif
diff --git a/utils/mf5to6/src/Preproc/ObsBlock.f90 b/utils/mf5to6/src/Preproc/ObsBlock.f90
index 39e6b64c984..59ca78fdff7 100644
--- a/utils/mf5to6/src/Preproc/ObsBlock.f90
+++ b/utils/mf5to6/src/Preproc/ObsBlock.f90
@@ -13,6 +13,7 @@ module ObsBlockModule
                                        GetObserveFromList, ConstructObservation
   use SimPHMFModule,             only: store_error, store_error_unit, ustop
   use UtilitiesModule,           only: CalcContribFactors
+  use GeomUtilModule,            only: get_node
   
   implicit none
   
@@ -74,7 +75,7 @@ subroutine process_block(this, insertLine, WriteBeginEnd, parser)
     type(BlockParserType), intent(inout) :: parser
     ! local
     integer :: ierr, ioutmf, layer
-    integer :: irow, jcol, lloc2, iadjrow, jadjcol, ncol, nrow
+    integer :: irow, jcol, lloc2, iadjrow, jadjcol, ncol, nrow, node
     double precision :: time, xcoord, xnode, ycoord, ynode, xoff, yoff
     double precision :: gridX, gridY
     character(len=MAXCHARLEN) :: ermsg
@@ -168,6 +169,9 @@ subroutine process_block(this, insertLine, WriteBeginEnd, parser)
         newobs%gridX = gridX
         newobs%gridY = gridY
         !
+        ! Determine node number based on indices
+        node = get_node(layer, irow, jcol, dis3d%nlay, dis3d%nrow, dis3d%ncol)
+        !
         ! Get coordinates of node and calculate X and Y offsets from node.
         call dis3d%get_node_coords_idx2(irow, jcol, xnode, ynode)
         xoff = gridX - xnode
@@ -187,14 +191,14 @@ subroutine process_block(this, insertLine, WriteBeginEnd, parser)
         if (.not. is_close(xoff, DZERO)) then
           if (xoff > DZERO) then
             if (jcol < ncol) then
-              if (dis3d%idomain(jcol+1, irow, layer) == 1) then
+              if (dis3d%idomain(node) == 1) then
                 jadjcol = jcol + 1
                 newobs%jcoladj = jadjcol
               endif
             endif
           else
             if (jcol > 1) then
-              if (dis3d%idomain(jcol-1, irow, layer) == 1) then
+              if (dis3d%idomain(node) == 1) then
                 jadjcol = jcol - 1
                 newobs%jcoladj = jadjcol
               endif
@@ -205,14 +209,14 @@ subroutine process_block(this, insertLine, WriteBeginEnd, parser)
         if (.not. is_close(yoff, DZERO)) then
           if (yoff > DZERO) then
             if (irow > 1) then
-              if (dis3d%idomain(jcol, irow-1, layer) == 1) then
+              if (dis3d%idomain(node) == 1) then
                 iadjrow = irow - 1
                 newobs%irowadj = iadjrow
               endif
             endif
           else
             if (irow < nrow) then
-              if (dis3d%idomain(jcol, irow+1, layer) == 1) then
+              if (dis3d%idomain(node) == 1) then
                 iadjrow = irow + 1
                 newobs%irowadj = iadjrow
               endif
@@ -223,7 +227,8 @@ subroutine process_block(this, insertLine, WriteBeginEnd, parser)
         ! If both jadjcol and iadjrow are not 0, and diagonal cell is active,
         ! assign newobs%idiag flag.
         if (jadjcol /= 0 .and. iadjrow /= 0) then
-          if (dis3d%idomain(jadjcol,iadjrow,layer) == 1) then
+          node = get_node(layer, iadjrow, jadjcol, dis3d%nlay, dis3d%nrow, dis3d%ncol)
+          if (dis3d%idomain(node) == 1) then
             newobs%idiag = 1
           endif
         endif