diff --git a/mpp/include/group_update_pack.inc b/mpp/include/group_update_pack.inc
index de08b89e56..98b7645836 100644
--- a/mpp/include/group_update_pack.inc
+++ b/mpp/include/group_update_pack.inc
@@ -17,9 +17,8 @@
 !***********************************************************************
 
 if( group%k_loop_inside ) then
-!$OMP parallel do default(none) shared(npack,group,ptr,nvector,ksize,buffer_start_pos) &
-!$OMP                          private(buffer_pos,pos,m,is,ie,js,je,rotation, &
-!$OMP                                  ptr_field, ptr_fieldx, ptr_fieldy,n,k)
+!$OMP parallel do default(none) shared(buffer,npack,group,nvector,ksize,buffer_start_pos) &
+!$OMP                          private(buffer_pos,pos,m,is,ie,js,je,rotation,n,k)
   do n = 1, npack
      buffer_pos = group%pack_buffer_pos(n) + buffer_start_pos
      pos  = buffer_pos
@@ -30,48 +29,48 @@ if( group%k_loop_inside ) then
         select case( rotation )
         case(ZERO)
            do l=1, group%nscalar  ! loop over number of fields
-              ptr_field = group%addrs_s(l)
               do k = 1, ksize
                  do j = js, je
                     do i = is, ie
                        pos = pos + 1
-                       buffer(pos) = field(i,j,k)
+                       !buffer(pos) = field(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_S, l, i, j, k)
                     end do
                  end do
               enddo
            enddo
         case( MINUS_NINETY )
            do l=1,group%nscalar  ! loop over number of fields
-              ptr_field = group%addrs_s(l)
               do k = 1, ksize
                  do i = is, ie
                     do j = je, js, -1
                        pos = pos + 1
-                       buffer(pos) = field(i,j,k)
+                       !buffer(pos) = field(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_S, l, i, j, k)
                     end do
                  end do
               end do
            end do
         case( NINETY )
            do l=1,group%nscalar  ! loop over number of fields
-              ptr_field = group%addrs_s(l)
               do k = 1, ksize
                  do i = ie, is, -1
                     do j = js, je
                        pos = pos + 1
-                       buffer(pos) = field(i,j,k)
+                       !buffer(pos) = field(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_S, l, i, j, k)
                     end do
                  end do
               end do
            end do
         case( ONE_HUNDRED_EIGHTY )
            do l=1,group%nscalar  ! loop over number of fields
-              ptr_field = group%addrs_s(l)
               do k = 1, ksize
                  do j = je, js, -1
                     do i = ie, is, -1
                        pos = pos + 1
-                       buffer(pos) = field(i,j,k)
+                       !buffer(pos) = field(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_S, l, i, j, k)
                     end do
                  end do
               end do
@@ -81,12 +80,12 @@ if( group%k_loop_inside ) then
         select case( rotation )
         case(ZERO)
            do l=1, nvector  ! loop over number of fields
-              ptr_fieldx = group%addrs_x(l)
               do k = 1, ksize
                  do j = js, je
                     do i = is, ie
                        pos = pos + 1
-                       buffer(pos) = fieldx(i,j,k)
+                       !buffer(pos) = fieldx(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_X, l, i, j, k)
                     end do
                  end do
               end do
@@ -94,24 +93,24 @@ if( group%k_loop_inside ) then
         case( MINUS_NINETY )
            if( BTEST(group%flags_v,SCALAR_BIT) ) then
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldy = group%addrs_y(l)
                  do k = 1, ksize
                     do i = is, ie
                        do j = je, js, -1
                           pos = pos + 1
-                          buffer(pos) = fieldy(i,j,k)
+                          !buffer(pos) = fieldy(i,j,k)
+                          buffer(pos) = GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                        end do
                     end do
                  end do
               end do
            else
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldy = group%addrs_y(l)
                  do k = 1, ksize
                     do i = is, ie
                        do j = je, js, -1
                           pos = pos + 1
-                          buffer(pos) = -fieldy(i,j,k)
+                          !buffer(pos) = -fieldy(i,j,k)
+                          buffer(pos) = - GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                        end do
                     end do
                  end do
@@ -119,12 +118,12 @@ if( group%k_loop_inside ) then
            end if
         case( NINETY )
            do l=1, nvector  ! loop over number of fields
-              ptr_fieldy = group%addrs_y(l)
               do k = 1, ksize
                  do i = ie, is, -1
                     do j = js, je
                        pos = pos + 1
-                       buffer(pos) = fieldy(i,j,k)
+                       !buffer(pos) = fieldy(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                     end do
                  end do
               end do
@@ -132,24 +131,24 @@ if( group%k_loop_inside ) then
         case( ONE_HUNDRED_EIGHTY )
            if( BTEST(group%flags_v,SCALAR_BIT) ) then
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldx = group%addrs_x(l)
                  do k = 1, ksize
                     do j = je, js, -1
                        do i = ie, is, -1
                           pos = pos + 1
-                          buffer(pos) =  fieldx(i,j,k)
+                          !buffer(pos) =  fieldx(i,j,k)
+                          buffer(pos) = GET_VALUE_ (group, FIELD_X, l, i, j, k)
                        end do
                     end do
                  end do
               end do
            else
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldx = group%addrs_x(l)
                  do k = 1, ksize
                     do j = je, js, -1
                        do i = ie, is, -1
                           pos = pos + 1
-                          buffer(pos) =  -fieldx(i,j,k)
+                          !buffer(pos) =  -fieldx(i,j,k)
+                          buffer(pos) = - GET_VALUE_ (group, FIELD_X, l, i, j, k)
                        end do
                     end do
                  end do
@@ -160,24 +159,24 @@ if( group%k_loop_inside ) then
         select case( rotation )
         case(ZERO)
            do l=1, nvector  ! loop over number of fields
-              ptr_fieldy = group%addrs_y(l)
               do k = 1, ksize
                  do j = js, je
                     do i = is, ie
                        pos = pos + 1
-                       buffer(pos) = fieldy(i,j,k)
+                       !buffer(pos) = fieldy(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                     end do
                  end do
               end do
            end do
         case( MINUS_NINETY )
            do l=1,nvector  ! loop over number of fields
-              ptr_fieldx = group%addrs_x(l)
               do k = 1, ksize
                  do i = is, ie
                     do j = je, js, -1
                        pos = pos + 1
-                       buffer(pos) = fieldx(i,j,k)
+                       !buffer(pos) = fieldx(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_X, l, i, j, k)
                     end do
                  end do
               end do
@@ -185,24 +184,24 @@ if( group%k_loop_inside ) then
         case( NINETY )
            if( BTEST(group%flags_v,SCALAR_BIT) ) then
               do l=1, nvector  ! loop over number of fields
-                 ptr_fieldx = group%addrs_x(l)
                  do k = 1, ksize
                     do i = ie, is, -1
                        do j = js, je
                           pos = pos + 1
-                          buffer(pos) = fieldx(i,j,k)
+                          !buffer(pos) = fieldx(i,j,k)
+                          buffer(pos) = GET_VALUE_ (group, FIELD_X, l, i, j, k)
                        end do
                     end do
                  end do
               end do
            else
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldx = group%addrs_x(l)
                  do k = 1, ksize
                     do i = ie, is, -1
                        do j = js, je
                           pos = pos + 1
-                          buffer(pos) = -fieldx(i,j,k)
+                          !buffer(pos) = -fieldx(i,j,k)
+                          buffer(pos) = - GET_VALUE_ (group, FIELD_X, l, i, j, k)
                        end do
                     end do
                  end do
@@ -211,24 +210,24 @@ if( group%k_loop_inside ) then
         case( ONE_HUNDRED_EIGHTY )
            if( BTEST(group%flags_v,SCALAR_BIT) ) then
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldy = group%addrs_y(l)
                  do k = 1, ksize
                     do j = je, js, -1
                        do i = ie, is, -1
                           pos = pos + 1
-                          buffer(pos) =  fieldy(i,j,k)
+                          !buffer(pos) =  fieldy(i,j,k)
+                          buffer(pos) = GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                        end do
                     end do
                  end do
               end do
            else
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldy = group%addrs_y(l)
                  do k = 1, ksize
                     do j = je, js, -1
                        do i = ie, is, -1
                           pos = pos + 1
-                          buffer(pos) =  -fieldy(i,j,k)
+                          !buffer(pos) =  -fieldy(i,j,k)
+                          buffer(pos) = - GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                        end do
                     end do
                  end do
@@ -238,9 +237,8 @@ if( group%k_loop_inside ) then
      endif
   enddo
 else
-!$OMP parallel do default(none) shared(npack,group,ptr,nvector,ksize,buffer_start_pos) &
-!$OMP                          private(buffer_pos,pos,m,is,ie,js,je,rotation, &
-!$OMP                                  ptr_field, ptr_fieldx, ptr_fieldy,n,k)
+!$OMP parallel do default(none) shared(buffer,npack,group,nvector,ksize,buffer_start_pos) &
+!$OMP                          private(buffer_pos,pos,m,is,ie,js,je,rotation,n,k)
   do nk = 1, npack*ksize
      n = (nk-1)/ksize + 1
      k = mod((nk-1), ksize) + 1
@@ -253,41 +251,41 @@ else
         select case( rotation )
         case(ZERO)
            do l=1, group%nscalar  ! loop over number of fields
-              ptr_field = group%addrs_s(l)
               do j = js, je
                  do i = is, ie
                     pos = pos + 1
-                    buffer(pos) = field(i,j,k)
+                    !buffer(pos) = field(i,j,k)
+                    buffer(pos) = GET_VALUE_ (group, FIELD_S, l, i, j, k)
                  end do
               end do
            enddo
         case( MINUS_NINETY )
            do l=1,group%nscalar  ! loop over number of fields
-              ptr_field = group%addrs_s(l)
               do i = is, ie
                  do j = je, js, -1
                     pos = pos + 1
-                    buffer(pos) = field(i,j,k)
+                    !buffer(pos) = field(i,j,k)
+                    buffer(pos) = GET_VALUE_ (group, FIELD_S, l, i, j, k)
                  end do
               end do
            end do
         case( NINETY )
            do l=1,group%nscalar  ! loop over number of fields
-              ptr_field = group%addrs_s(l)
               do i = ie, is, -1
                  do j = js, je
                     pos = pos + 1
-                    buffer(pos) = field(i,j,k)
+                    !buffer(pos) = field(i,j,k)
+                    buffer(pos) = GET_VALUE_ (group, FIELD_S, l, i, j, k)
                  end do
               end do
            end do
         case( ONE_HUNDRED_EIGHTY )
            do l=1,group%nscalar  ! loop over number of fields
-              ptr_field = group%addrs_s(l)
               do j = je, js, -1
                  do i = ie, is, -1
                     pos = pos + 1
-                    buffer(pos) = field(i,j,k)
+                    !buffer(pos) = field(i,j,k)
+                    buffer(pos) = GET_VALUE_ (group, FIELD_S, l, i, j, k)
                  end do
               end do
            end do
@@ -296,64 +294,64 @@ else
         select case( rotation )
         case(ZERO)
            do l=1, nvector  ! loop over number of fields
-              ptr_fieldx = group%addrs_x(l)
               do j = js, je
                  do i = is, ie
                     pos = pos + 1
-                    buffer(pos) = fieldx(i,j,k)
+                    !buffer(pos) = fieldx(i,j,k)
+                    buffer(pos) = GET_VALUE_ (group, FIELD_X, l, i, j, k)
                  end do
               end do
            end do
         case( MINUS_NINETY )
            if( BTEST(group%flags_v,SCALAR_BIT) ) then
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldy = group%addrs_y(l)
                  do i = is, ie
                     do j = je, js, -1
                        pos = pos + 1
-                       buffer(pos) = fieldy(i,j,k)
+                       !buffer(pos) = fieldy(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                     end do
                  end do
               end do
            else
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldy = group%addrs_y(l)
                  do i = is, ie
                     do j = je, js, -1
                        pos = pos + 1
-                       buffer(pos) = -fieldy(i,j,k)
+                       !buffer(pos) = -fieldy(i,j,k)
+                       buffer(pos) = - GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                     end do
                  end do
               end do
            end if
         case( NINETY )
            do l=1, nvector  ! loop over number of fields
-              ptr_fieldy = group%addrs_y(l)
               do i = ie, is, -1
                  do j = js, je
                     pos = pos + 1
-                    buffer(pos) = fieldy(i,j,k)
+                    !buffer(pos) = fieldy(i,j,k)
+                    buffer(pos) = GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                  end do
               end do
            end do
         case( ONE_HUNDRED_EIGHTY )
            if( BTEST(group%flags_v,SCALAR_BIT) ) then
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldx = group%addrs_x(l)
                  do j = je, js, -1
                     do i = ie, is, -1
                        pos = pos + 1
-                       buffer(pos) =  fieldx(i,j,k)
+                       !buffer(pos) =  fieldx(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_X, l, i, j, k)
                     end do
                  end do
               end do
            else
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldx = group%addrs_x(l)
                  do j = je, js, -1
                     do i = ie, is, -1
                        pos = pos + 1
-                       buffer(pos) =  -fieldx(i,j,k)
+                       !buffer(pos) =  -fieldx(i,j,k)
+                       buffer(pos) = - GET_VALUE_ (group, FIELD_X, l, i, j, k)
                     end do
                  end do
               end do
@@ -363,42 +361,42 @@ else
         select case( rotation )
         case(ZERO)
            do l=1, nvector  ! loop over number of fields
-              ptr_fieldy = group%addrs_y(l)
               do j = js, je
                  do i = is, ie
                     pos = pos + 1
-                    buffer(pos) = fieldy(i,j,k)
+                    !buffer(pos) = fieldy(i,j,k)
+                    buffer(pos) = GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                  end do
               end do
            end do
         case( MINUS_NINETY )
            do l=1,nvector  ! loop over number of fields
-              ptr_fieldx = group%addrs_x(l)
               do i = is, ie
                  do j = je, js, -1
                     pos = pos + 1
-                    buffer(pos) = fieldx(i,j,k)
+                    !buffer(pos) = fieldx(i,j,k)
+                    buffer(pos) = GET_VALUE_ (group, FIELD_X, l, i, j, k)
                  end do
               end do
            end do
         case( NINETY )
            if( BTEST(group%flags_v,SCALAR_BIT) ) then
               do l=1, nvector  ! loop over number of fields
-                 ptr_fieldx = group%addrs_x(l)
                  do i = ie, is, -1
                     do j = js, je
                        pos = pos + 1
-                       buffer(pos) = fieldx(i,j,k)
+                       !buffer(pos) = fieldx(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_X, l, i, j, k)
                     end do
                  end do
               end do
            else
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldx = group%addrs_x(l)
                  do i = ie, is, -1
                     do j = js, je
                        pos = pos + 1
-                       buffer(pos) = -fieldx(i,j,k)
+                       !buffer(pos) = -fieldx(i,j,k)
+                       buffer(pos) = - GET_VALUE_ (group, FIELD_X, l, i, j, k)
                     end do
                  end do
               end do
@@ -406,21 +404,21 @@ else
         case( ONE_HUNDRED_EIGHTY )
            if( BTEST(group%flags_v,SCALAR_BIT) ) then
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldy = group%addrs_y(l)
                  do j = je, js, -1
                     do i = ie, is, -1
                        pos = pos + 1
-                       buffer(pos) =  fieldy(i,j,k)
+                       !buffer(pos) =  fieldy(i,j,k)
+                       buffer(pos) = GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                     end do
                  end do
               end do
            else
               do l=1,nvector  ! loop over number of fields
-                 ptr_fieldy = group%addrs_y(l)
                  do j = je, js, -1
                     do i = ie, is, -1
                        pos = pos + 1
-                       buffer(pos) =  -fieldy(i,j,k)
+                       !buffer(pos) =  -fieldy(i,j,k)
+                       buffer(pos) = - GET_VALUE_ (group, FIELD_Y, l, i, j, k)
                     end do
                  end do
               end do
diff --git a/mpp/include/group_update_unpack.inc b/mpp/include/group_update_unpack.inc
index 49fb2555ce..06c8a36bed 100644
--- a/mpp/include/group_update_unpack.inc
+++ b/mpp/include/group_update_unpack.inc
@@ -17,9 +17,8 @@
 !***********************************************************************
 
 if( group%k_loop_inside ) then
-!$OMP parallel do default(none) shared(nunpack,group,nscalar,ptr,nvector,ksize,buffer_start_pos) &
-!$OMP                          private(buffer_pos,pos,m,is, ie, js, je,rotation, &
-!$OMP                                  ptr_field, ptr_fieldx, ptr_fieldy, n,k )
+!$OMP parallel do default(none) shared(buffer,nunpack,group,nscalar,nvector,ksize,buffer_start_pos) &
+!$OMP                          private(buffer_pos,pos,m,is, ie, js, je,rotation,n,k )
   do n = nunpack, 1, -1
      buffer_pos = group%unpack_buffer_pos(n) + buffer_start_pos
      pos = buffer_pos
@@ -27,36 +26,36 @@ if( group%k_loop_inside ) then
      js = group%unpack_js(n); je = group%unpack_je(n)
      if( group%unpack_type(n) == FIELD_S ) then
         do l=1,nscalar  ! loop over number of fields
-           ptr_field = group%addrs_s(l)
            do k = 1, ksize
               do j = js, je
                  do i = is, ie
                     pos = pos + 1
-                    field(i,j,k) = buffer(pos)
+                    !field(i,j,k) = buffer(pos)
+                    call SET_VALUE_ (group, FIELD_S, l, i, j, k, buffer(pos))
                  end do
               end do
            end do
         end do
      else if( group%unpack_type(n) == FIELD_X ) then
         do l=1,nvector  ! loop over number of fields
-           ptr_fieldx = group%addrs_x(l)
            do k = 1, ksize
               do j = js, je
                  do i = is, ie
                     pos = pos + 1
-                    fieldx(i,j,k) = buffer(pos)
+                    !fieldx(i,j,k) = buffer(pos)
+                    call SET_VALUE_ (group, FIELD_X, l, i, j, k, buffer(pos))
                  end do
               end do
            end do
         end do
      else if( group%unpack_type(n) == FIELD_Y ) then
         do l=1,nvector  ! loop over number of fields
-           ptr_fieldy = group%addrs_y(l)
            do k = 1, ksize
               do j = js, je
                  do i = is, ie
                     pos = pos + 1
-                    fieldy(i,j,k) = buffer(pos)
+                    !fieldy(i,j,k) = buffer(pos)
+                    call SET_VALUE_ (group, FIELD_Y, l, i, j, k, buffer(pos))
                  end do
               end do
            end do
@@ -64,9 +63,8 @@ if( group%k_loop_inside ) then
      endif
   enddo
 else
-!$OMP parallel do default(none) shared(nunpack,group,nscalar,ptr,nvector,ksize,buffer_start_pos) &
-!$OMP                          private(buffer_pos,pos,m,is, ie, js, je,rotation, &
-!$OMP                                  ptr_field, ptr_fieldx, ptr_fieldy,n,k)
+!$OMP parallel do default(none) shared(buffer,nunpack,group,nscalar,nvector,ksize,buffer_start_pos) &
+!$OMP                          private(buffer_pos,pos,m,is, ie, js, je,rotation,n,k)
   do nk = nunpack*ksize, 1, -1
      n = (nk-1)/ksize + 1
      k = mod((nk-1), ksize) + 1
@@ -76,31 +74,31 @@ else
      js = group%unpack_js(n); je = group%unpack_je(n)
      if( group%unpack_type(n) == FIELD_S ) then
         do l=1,nscalar  ! loop over number of fields
-           ptr_field = group%addrs_s(l)
            do j = js, je
               do i = is, ie
                  pos = pos + 1
-                 field(i,j,k) = buffer(pos)
+                 !field(i,j,k) = buffer(pos)
+                 call SET_VALUE_ (group, FIELD_S, l, i, j, k, buffer(pos))
               end do
            end do
         end do
      else if( group%unpack_type(n) == FIELD_X ) then
         do l=1,nvector  ! loop over number of fields
-           ptr_fieldx = group%addrs_x(l)
            do j = js, je
               do i = is, ie
                  pos = pos + 1
-                 fieldx(i,j,k) = buffer(pos)
+                 !fieldx(i,j,k) = buffer(pos)
+                 call SET_VALUE_ (group, FIELD_X, l, i, j, k, buffer(pos))
               end do
            end do
         end do
      else if( group%unpack_type(n) == FIELD_Y ) then
         do l=1,nvector  ! loop over number of fields
-           ptr_fieldy = group%addrs_y(l)
            do j = js, je
               do i = is, ie
                  pos = pos + 1
-                 fieldy(i,j,k) = buffer(pos)
+                 !fieldy(i,j,k) = buffer(pos)
+                 call SET_VALUE_ (group, FIELD_Y, l, i, j, k, buffer(pos))
               end do
            end do
         end do
diff --git a/mpp/include/mpp_data_mpi.inc b/mpp/include/mpp_data_mpi.inc
index b717c78e2b..466a533f8d 100644
--- a/mpp/include/mpp_data_mpi.inc
+++ b/mpp/include/mpp_data_mpi.inc
@@ -38,8 +38,8 @@ integer, parameter :: mpp_from_pe = -999, ptr_from = -999
 !-------------------------------------------------------------------!
 ! The following data is used in mpp_domains_mod and its components  !
 !-------------------------------------------------------------------!
-real(r8_kind), allocatable :: mpp_domains_stack(:) !< stack used to hold data for domain operations
-real(r8_kind), allocatable :: mpp_domains_stack_nonblock(:) !< stack used for non-blocking domain operations
+real(r8_kind), allocatable, target :: mpp_domains_stack(:) !< stack used to hold data for domain operations
+real(r8_kind), allocatable, target :: mpp_domains_stack_nonblock(:) !< stack used for non-blocking domain operations
 !--- some dummy variables with dummy values that will never be used
 integer, parameter :: ptr_domains_stack = -999
 integer, parameter :: ptr_domains_stack_nonblock = -999
diff --git a/mpp/include/mpp_data_nocomm.inc b/mpp/include/mpp_data_nocomm.inc
index f7de0670e4..62d4bd8065 100644
--- a/mpp/include/mpp_data_nocomm.inc
+++ b/mpp/include/mpp_data_nocomm.inc
@@ -38,8 +38,8 @@ integer, parameter :: mpp_from_pe = -999, ptr_from = -999
 !-------------------------------------------------------------------!
 ! The following data is used in mpp_domains_mod and its components  !
 !-------------------------------------------------------------------!
-real(r8_kind), allocatable :: mpp_domains_stack(:)
-real(r8_kind), allocatable :: mpp_domains_stack_nonblock(:)
+real(r8_kind), allocatable, target :: mpp_domains_stack(:)
+real(r8_kind), allocatable, target :: mpp_domains_stack_nonblock(:)
 !--- some dummy variables with dummy values that will never be used
 integer, parameter :: ptr_domains_stack = -999
 integer, parameter :: ptr_domains_stack_nonblock = -999
diff --git a/mpp/include/mpp_domains_misc.inc b/mpp/include/mpp_domains_misc.inc
index 82d60358db..6bff08a870 100644
--- a/mpp/include/mpp_domains_misc.inc
+++ b/mpp/include/mpp_domains_misc.inc
@@ -1956,70 +1956,46 @@ end subroutine init_nonblock_type
 #define MPP_TYPE_ real(r8_kind)
 #undef MPI_TYPE_
 #define MPI_TYPE_ MPI_REAL8
-#undef MPP_CREATE_GROUP_UPDATE_2D_
-#define MPP_CREATE_GROUP_UPDATE_2D_ mpp_create_group_update_r8_2d
-#undef MPP_CREATE_GROUP_UPDATE_3D_
-#define MPP_CREATE_GROUP_UPDATE_3D_ mpp_create_group_update_r8_3d
-#undef MPP_CREATE_GROUP_UPDATE_4D_
-#define MPP_CREATE_GROUP_UPDATE_4D_ mpp_create_group_update_r8_4d
-#undef MPP_CREATE_GROUP_UPDATE_2D_V_
-#define MPP_CREATE_GROUP_UPDATE_2D_V_ mpp_create_group_update_r8_2dv
-#undef MPP_CREATE_GROUP_UPDATE_3D_V_
-#define MPP_CREATE_GROUP_UPDATE_3D_V_ mpp_create_group_update_r8_3dv
-#undef MPP_CREATE_GROUP_UPDATE_4D_V_
-#define MPP_CREATE_GROUP_UPDATE_4D_V_ mpp_create_group_update_r8_4dv
+#undef MPP_CREATE_GROUP_UPDATE_
+#define MPP_CREATE_GROUP_UPDATE_ mpp_create_group_update_r8
+#undef MPP_CREATE_GROUP_UPDATE_V_
+#define MPP_CREATE_GROUP_UPDATE_V_ mpp_create_group_update_r8_v
 #undef MPP_DO_GROUP_UPDATE_
 #define MPP_DO_GROUP_UPDATE_ mpp_do_group_update_r8
 #undef MPP_START_GROUP_UPDATE_
 #define MPP_START_GROUP_UPDATE_ mpp_start_group_update_r8
 #undef MPP_COMPLETE_GROUP_UPDATE_
 #define MPP_COMPLETE_GROUP_UPDATE_ mpp_complete_group_update_r8
-#undef MPP_RESET_GROUP_UPDATE_FIELD_2D_
-#define MPP_RESET_GROUP_UPDATE_FIELD_2D_ mpp_reset_group_update_field_r8_2d
-#undef MPP_RESET_GROUP_UPDATE_FIELD_3D_
-#define MPP_RESET_GROUP_UPDATE_FIELD_3D_ mpp_reset_group_update_field_r8_3d
-#undef MPP_RESET_GROUP_UPDATE_FIELD_4D_
-#define MPP_RESET_GROUP_UPDATE_FIELD_4D_ mpp_reset_group_update_field_r8_4d
-#undef MPP_RESET_GROUP_UPDATE_FIELD_2D_V_
-#define MPP_RESET_GROUP_UPDATE_FIELD_2D_V_ mpp_reset_group_update_field_r8_2dv
-#undef MPP_RESET_GROUP_UPDATE_FIELD_3D_V_
-#define MPP_RESET_GROUP_UPDATE_FIELD_3D_V_ mpp_reset_group_update_field_r8_3dv
-#undef MPP_RESET_GROUP_UPDATE_FIELD_4D_V_
-#define MPP_RESET_GROUP_UPDATE_FIELD_4D_V_ mpp_reset_group_update_field_r8_4dv
+#undef MPP_RESET_GROUP_UPDATE_FIELD_
+#define MPP_RESET_GROUP_UPDATE_FIELD_ mpp_reset_group_update_field_r8
+#undef MPP_RESET_GROUP_UPDATE_FIELD_V_
+#define MPP_RESET_GROUP_UPDATE_FIELD_V_ mpp_reset_group_update_field_r8_v
+#undef GET_VALUE_
+#define GET_VALUE_ get_value_r8
+#undef SET_VALUE_
+#define SET_VALUE_ set_value_r8
 #include <mpp_group_update.fh>
 
 #undef MPP_TYPE_
 #define MPP_TYPE_ real(r4_kind)
 #undef MPI_TYPE_
 #define MPI_TYPE_ MPI_REAL4
-#undef MPP_CREATE_GROUP_UPDATE_2D_
-#define MPP_CREATE_GROUP_UPDATE_2D_ mpp_create_group_update_r4_2d
-#undef MPP_CREATE_GROUP_UPDATE_3D_
-#define MPP_CREATE_GROUP_UPDATE_3D_ mpp_create_group_update_r4_3d
-#undef MPP_CREATE_GROUP_UPDATE_4D_
-#define MPP_CREATE_GROUP_UPDATE_4D_ mpp_create_group_update_r4_4d
-#undef MPP_CREATE_GROUP_UPDATE_2D_V_
-#define MPP_CREATE_GROUP_UPDATE_2D_V_ mpp_create_group_update_r4_2dv
-#undef MPP_CREATE_GROUP_UPDATE_3D_V_
-#define MPP_CREATE_GROUP_UPDATE_3D_V_ mpp_create_group_update_r4_3dv
-#undef MPP_CREATE_GROUP_UPDATE_4D_V_
-#define MPP_CREATE_GROUP_UPDATE_4D_V_ mpp_create_group_update_r4_4dv
+#undef MPP_CREATE_GROUP_UPDATE_
+#define MPP_CREATE_GROUP_UPDATE_ mpp_create_group_update_r4
+#undef MPP_CREATE_GROUP_UPDATE_V_
+#define MPP_CREATE_GROUP_UPDATE_V_ mpp_create_group_update_r4_v
 #undef MPP_DO_GROUP_UPDATE_
 #define MPP_DO_GROUP_UPDATE_ mpp_do_group_update_r4
 #undef MPP_START_GROUP_UPDATE_
 #define MPP_START_GROUP_UPDATE_ mpp_start_group_update_r4
 #undef MPP_COMPLETE_GROUP_UPDATE_
 #define MPP_COMPLETE_GROUP_UPDATE_ mpp_complete_group_update_r4
-#undef MPP_RESET_GROUP_UPDATE_FIELD_2D_
-#define MPP_RESET_GROUP_UPDATE_FIELD_2D_ mpp_reset_group_update_field_r4_2d
-#undef MPP_RESET_GROUP_UPDATE_FIELD_3D_
-#define MPP_RESET_GROUP_UPDATE_FIELD_3D_ mpp_reset_group_update_field_r4_3d
-#undef MPP_RESET_GROUP_UPDATE_FIELD_4D_
-#define MPP_RESET_GROUP_UPDATE_FIELD_4D_ mpp_reset_group_update_field_r4_4d
-#undef MPP_RESET_GROUP_UPDATE_FIELD_2D_V_
-#define MPP_RESET_GROUP_UPDATE_FIELD_2D_V_ mpp_reset_group_update_field_r4_2dv
-#undef MPP_RESET_GROUP_UPDATE_FIELD_3D_V_
-#define MPP_RESET_GROUP_UPDATE_FIELD_3D_V_ mpp_reset_group_update_field_r4_3dv
-#undef MPP_RESET_GROUP_UPDATE_FIELD_4D_V_
-#define MPP_RESET_GROUP_UPDATE_FIELD_4D_V_ mpp_reset_group_update_field_r4_4dv
+#undef MPP_RESET_GROUP_UPDATE_FIELD_
+#define MPP_RESET_GROUP_UPDATE_FIELD_ mpp_reset_group_update_field_r4
+#undef MPP_RESET_GROUP_UPDATE_FIELD_V_
+#define MPP_RESET_GROUP_UPDATE_FIELD_V_ mpp_reset_group_update_field_r4_v
+#undef GET_VALUE_
+#define GET_VALUE_ get_value_r4
+#undef SET_VALUE_
+#define SET_VALUE_ set_value_r4
 #include <mpp_group_update.fh>
diff --git a/mpp/include/mpp_group_update.fh b/mpp/include/mpp_group_update.fh
index 8373ab0290..c84325ac6c 100644
--- a/mpp/include/mpp_group_update.fh
+++ b/mpp/include/mpp_group_update.fh
@@ -18,32 +18,15 @@
 ! -*-f90-*-
 !> @addtogroup mpp_domains_mod
 !> @{
-subroutine MPP_CREATE_GROUP_UPDATE_2D_(group, field, domain, flags, position, &
-     whalo, ehalo, shalo, nhalo)
+subroutine MPP_CREATE_GROUP_UPDATE_(group, field, domain, flags, position, whalo, ehalo, shalo, nhalo, xdim, ydim)
   type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,            intent(inout)        :: field(:,:)
-  type(domain2D),       intent(inout)        :: domain
-  integer,              intent(in), optional :: flags
-  integer,              intent(in), optional :: position
-  integer,              intent(in), optional :: whalo, ehalo, shalo, nhalo
-
-  MPP_TYPE_ :: field3D(size(field,1),size(field,2),1)
-  pointer( ptr, field3D )
-  ptr = LOC(field)
-
-  call mpp_create_group_update(group, field3D, domain, flags, position, whalo, ehalo, shalo, nhalo)
-
-  return
-
-end subroutine MPP_CREATE_GROUP_UPDATE_2D_
-
-subroutine MPP_CREATE_GROUP_UPDATE_3D_(group, field, domain, flags, position, whalo, ehalo, shalo, nhalo)
-  type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,        intent(inout)        :: field(:,:,:)
+  MPP_TYPE_,        intent(inout), target :: field(..)
   type(domain2D),   intent(inout)        :: domain
   integer,          intent(in), optional :: flags
   integer,          intent(in), optional :: position
   integer,          intent(in), optional :: whalo, ehalo, shalo, nhalo ! specify halo region to be updated.
+  integer,          intent(in), optional :: xdim, ydim !< Domain-decomposed dimensions. If not specified, it is assumed
+                                                       !! that the first two dimensions are x and y respectively.
 
   integer          :: update_position, update_whalo, update_ehalo, update_shalo, update_nhalo
   integer          :: update_flags, isize, jsize, ksize
@@ -51,11 +34,18 @@ subroutine MPP_CREATE_GROUP_UPDATE_3D_(group, field, domain, flags, position, wh
   character(len=3) :: text
   logical          :: set_mismatch, update_edge_only
   logical          :: recv(8)
+  integer          :: ix, iy
 
   if(group%initialized) then
       call mpp_error(FATAL, "MPP_CREATE_GROUP_UPDATE_3D: group is already initialized")
   endif
 
+  ix = 1
+  iy = 2
+
+  if (present(xdim)) ix=xdim
+  if (present(ydim)) iy=ydim
+
   if(present(whalo)) then
      update_whalo = whalo
      if(abs(update_whalo) > domain%whalo ) call mpp_error(FATAL, "MPP_CREATE_GROUP_UPDATE: "// &
@@ -107,9 +97,10 @@ subroutine MPP_CREATE_GROUP_UPDATE_3D_(group, field, domain, flags, position, wh
      call mpp_error(FATAL,'MPP_CREATE_GROUP_UPDATE: MAX_DOMAIN_FIELDS='//text//' exceeded for group update.' )
   endif
 
-  isize = size(field,1); jsize=size(field,2); ksize = size(field,3)
+  isize = size(field,ix); jsize=size(field,iy); ksize = size(field) / (isize*jsize)
+
+  group%addrs_s(nscalar) = c_loc(field)
 
-  group%addrs_s(nscalar) = LOC(field)
   if( group%nscalar == 1 ) then
      group%flags_s  = update_flags
      group%whalo_s  = update_whalo
@@ -120,8 +111,13 @@ subroutine MPP_CREATE_GROUP_UPDATE_3D_(group, field, domain, flags, position, wh
      group%isize_s  = isize
      group%jsize_s  = jsize
      group%ksize_s  = ksize
+     group%ix_s     = ix
+     group%iy_s     = iy
      call mpp_get_memory_domain(domain, group%is_s, group%ie_s, group%js_s, group%je_s, position=position)
 
+     allocate (group%shape_s(rank(field)))
+     group%shape_s = shape(field)
+
      update_edge_only = BTEST(update_flags, EDGEONLY)
      recv(1) = BTEST(update_flags,EAST)
      recv(3) = BTEST(update_flags,SOUTH)
@@ -163,61 +159,17 @@ subroutine MPP_CREATE_GROUP_UPDATE_3D_(group, field, domain, flags, position, wh
      endif
   endif
 
-  return
-
-end subroutine MPP_CREATE_GROUP_UPDATE_3D_
-
-
-subroutine MPP_CREATE_GROUP_UPDATE_4D_(group, field, domain, flags, position, &
-     whalo, ehalo, shalo, nhalo)
-  type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,            intent(inout)        :: field(:,:,:,:)
-  type(domain2D),       intent(inout)        :: domain
-  integer,              intent(in), optional :: flags
-  integer,              intent(in), optional :: position
-  integer,              intent(in), optional :: whalo, ehalo, shalo, nhalo
-
-  MPP_TYPE_ :: field3D(size(field,1),size(field,2),size(field,3)*size(field,4))
-  pointer( ptr, field3D )
-  ptr = LOC(field)
-
-  call mpp_create_group_update(group, field3D, domain, flags, position, whalo, ehalo, shalo, nhalo)
-
-  return
-
-end subroutine MPP_CREATE_GROUP_UPDATE_4D_
+end subroutine MPP_CREATE_GROUP_UPDATE_
 
-subroutine MPP_CREATE_GROUP_UPDATE_2D_V_( group, fieldx, fieldy, domain, flags, gridtype, &
-                                          whalo, ehalo, shalo, nhalo)
 
+subroutine MPP_CREATE_GROUP_UPDATE_V_( group, fieldx, fieldy, domain, flags, gridtype, &
+                                          whalo, ehalo, shalo, nhalo, xdim, ydim)
   type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,        intent(inout)            :: fieldx(:,:), fieldy(:,:)
-  type(domain2D),   intent(inout)            :: domain
-  integer,          intent(in),     optional :: flags, gridtype
-  integer,          intent(in),     optional :: whalo, ehalo, shalo, nhalo
-  MPP_TYPE_ :: field3Dx(size(fieldx,1),size(fieldx,2),1)
-  MPP_TYPE_ :: field3Dy(size(fieldy,1),size(fieldy,2),1)
-  pointer( ptrx, field3Dx )
-  pointer( ptry, field3Dy )
-  ptrx = LOC(fieldx)
-  ptry = LOC(fieldy)
-
-  call  mpp_create_group_update(group, field3Dx, field3Dy, domain, flags, gridtype, &
-                                          whalo, ehalo, shalo, nhalo)
-
-  return
-
-end subroutine MPP_CREATE_GROUP_UPDATE_2D_V_
-
-
-
-subroutine MPP_CREATE_GROUP_UPDATE_3D_V_( group, fieldx, fieldy, domain, flags, gridtype, &
-                                          whalo, ehalo, shalo, nhalo)
-  type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,        intent(inout)        :: fieldx(:,:,:), fieldy(:,:,:)
+  MPP_TYPE_,        intent(inout), target :: fieldx(..), fieldy(..)
   type(domain2D),   intent(inout)        :: domain
   integer,          intent(in), optional :: flags, gridtype
   integer,          intent(in), optional :: whalo, ehalo, shalo, nhalo
+  integer,          intent(in), optional :: xdim, ydim
 
   integer          :: update_whalo, update_ehalo, update_shalo, update_nhalo
   integer          :: update_flags, isize_x, jsize_x, ksize_x, isize_y, jsize_y, ksize_y
@@ -225,12 +177,19 @@ subroutine MPP_CREATE_GROUP_UPDATE_3D_V_( group, fieldx, fieldy, domain, flags,
   character(len=3) :: text
   logical          :: set_mismatch, update_edge_only
   logical          :: recv(8)
+  integer          :: ix, iy
 
 
   if(group%initialized) then
      call mpp_error(FATAL, "MPP_CREATE_GROUP_UPDATE_V: group is already initialized")
   endif
 
+  ix = 1
+  iy = 2
+
+  if (present(xdim)) ix=xdim
+  if (present(ydim)) iy=ydim
+
   if(present(whalo)) then
      update_whalo = whalo
      if(abs(update_whalo) > domain%whalo ) call mpp_error(FATAL, "MPP_CREATE_GROUP_UPDATE_V: "// &
@@ -284,14 +243,14 @@ subroutine MPP_CREATE_GROUP_UPDATE_3D_V_( group, fieldx, fieldy, domain, flags,
      call mpp_error(FATAL,'MPP_CREATE_GROUP_UPDATE_V: MAX_DOMAIN_FIELDS='//text//' exceeded for group update.' )
   endif
 
-  isize_x = size(fieldx,1); jsize_x = size(fieldx,2); ksize_x = size(fieldx,3)
-  isize_y = size(fieldy,1); jsize_y = size(fieldy,2); ksize_y = size(fieldy,3)
+  isize_x = size(fieldx,ix); jsize_x = size(fieldx,iy); ksize_x = size(fieldx) / (isize_x*jsize_x)
+  isize_y = size(fieldy,ix); jsize_y = size(fieldy,iy); ksize_y = size(fieldy) / (isize_y*jsize_y)
 
   if(ksize_x .NE. ksize_y) call mpp_error(FATAL,  &
            'MPP_CREATE_GROUP_UPDATE_V: mismatch of ksize between fieldx and fieldy')
 
-  group%addrs_x(nvector) = LOC(fieldx)
-  group%addrs_y(nvector) = LOC(fieldy)
+  group%addrs_x(nvector) = c_loc(fieldx)
+  group%addrs_y(nvector) = c_loc(fieldy)
 
   if( group%nvector == 1 ) then
      group%flags_v  = update_flags
@@ -305,6 +264,8 @@ subroutine MPP_CREATE_GROUP_UPDATE_3D_V_( group, fieldx, fieldy, domain, flags,
      group%isize_y  = isize_y
      group%jsize_y  = jsize_y
      group%ksize_v  = ksize_x
+     group%ix_v     = ix
+     group%iy_v     = iy
      update_edge_only = BTEST(update_flags, EDGEONLY)
      group%nonsym_edge = .false.
 
@@ -371,6 +332,12 @@ subroutine MPP_CREATE_GROUP_UPDATE_3D_V_( group, fieldx, fieldy, domain, flags,
 
      call mpp_get_memory_domain(domain, group%is_x, group%ie_x, group%js_x, group%je_x, position=position_x)
      call mpp_get_memory_domain(domain, group%is_y, group%ie_y, group%js_y, group%je_y, position=position_y)
+
+     allocate (group%shape_x(rank(fieldx)))
+     group%shape_x = shape(fieldx)
+
+     allocate (group%shape_y(rank(fieldy)))
+     group%shape_y = shape(fieldy)
   else
      set_mismatch = .false.
      set_mismatch = set_mismatch .OR. (group%flags_v  .NE. update_flags)
@@ -391,31 +358,7 @@ subroutine MPP_CREATE_GROUP_UPDATE_3D_V_( group, fieldx, fieldy, domain, flags,
      endif
   endif
 
-  return
-
-end subroutine MPP_CREATE_GROUP_UPDATE_3D_V_
-
-subroutine MPP_CREATE_GROUP_UPDATE_4D_V_( group, fieldx, fieldy, domain, flags, gridtype, &
-                                          whalo, ehalo, shalo, nhalo)
-
-  type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,        intent(inout)            :: fieldx(:,:,:,:), fieldy(:,:,:,:)
-  type(domain2D),   intent(inout)            :: domain
-  integer,          intent(in),     optional :: flags, gridtype
-  integer,          intent(in),     optional :: whalo, ehalo, shalo, nhalo
-  MPP_TYPE_ :: field3Dx(size(fieldx,1),size(fieldx,2),size(fieldx,3)*size(fieldx,4))
-  MPP_TYPE_ :: field3Dy(size(fieldy,1),size(fieldy,2),size(fieldy,3)*size(fieldy,4))
-  pointer( ptrx, field3Dx )
-  pointer( ptry, field3Dy )
-  ptrx = LOC(fieldx)
-  ptry = LOC(fieldy)
-
-  call  mpp_create_group_update(group, field3Dx, field3Dy, domain, flags, gridtype, &
-                                          whalo, ehalo, shalo, nhalo)
-
-  return
-
-end subroutine MPP_CREATE_GROUP_UPDATE_4D_V_
+end subroutine MPP_CREATE_GROUP_UPDATE_V_
 
 
 subroutine MPP_DO_GROUP_UPDATE_(group, domain, d_type)
@@ -423,7 +366,7 @@ subroutine MPP_DO_GROUP_UPDATE_(group, domain, d_type)
   type(domain2D),              intent(inout) :: domain
   MPP_TYPE_,                   intent(in)    :: d_type
 
-  integer   :: nscalar, nvector, nlist
+  integer   :: nscalar, nvector
   logical   :: recv_y(8)
   integer   :: nsend, nrecv, flags_v
   integer   :: msgsize
@@ -432,19 +375,13 @@ subroutine MPP_DO_GROUP_UPDATE_(group, domain, d_type)
   integer   :: n, l, m, i, j, k, buffer_start_pos, nk
   integer   :: shift, gridtype, midpoint
   integer   :: npack, nunpack, rotation, isd
+  type(c_ptr) :: stack_cptr !< Workaround for GFortran bug
 
-  MPP_TYPE_ :: buffer(mpp_domains_stack_size)
-  MPP_TYPE_ :: field (group%is_s:group%ie_s,group%js_s:group%je_s, group%ksize_s)
-  MPP_TYPE_ :: fieldx(group%is_x:group%ie_x,group%js_x:group%je_x, group%ksize_v)
-  MPP_TYPE_ :: fieldy(group%is_y:group%ie_y,group%js_y:group%je_y, group%ksize_v)
-  pointer(ptr, buffer )
-  pointer(ptr_field, field)
-  pointer(ptr_fieldx, fieldx)
-  pointer(ptr_fieldy, fieldy)
+  MPP_TYPE_, pointer :: buffer(:)
+  MPP_TYPE_, parameter :: zero_ = 0.
 
   nscalar = group%nscalar
   nvector = group%nvector
-  nlist   = size(domain%list(:))
   gridtype = group%gridtype
 
   !--- ksize_s must equal ksize_v
@@ -462,7 +399,8 @@ subroutine MPP_DO_GROUP_UPDATE_(group, domain, d_type)
   endif
   if(nvector > 0) recv_y = group%recv_y
 
-  ptr = LOC(mpp_domains_stack)
+  stack_cptr = c_loc(mpp_domains_stack(1))
+  call c_f_pointer(stack_cptr, buffer, [mpp_domains_stack_size])
 
   !--- set reset_index_s and reset_index_v to 0
   group%reset_index_s = 0
@@ -535,11 +473,9 @@ subroutine MPP_DO_GROUP_UPDATE_(group, domain, d_type)
            do i = is ,ie, midpoint
               if( domain%x(1)%domain_data%begin.LE.i .AND. i.LE. domain%x(1)%domain_data%end+shift )then
                  do l=1,nvector
-                    ptr_fieldx = group%addrs_x(l)
-                    ptr_fieldy = group%addrs_y(l)
                     do k = 1,ksize
-                       fieldx(i,j,k) = 0.
-                       fieldy(i,j,k) = 0.
+                       call SET_VALUE_ (group, FIELD_X, l, i, j, k, zero_)
+                       call SET_VALUE_ (group, FIELD_Y, l, i, j, k, zero_)
                     end do
                  end do
               end if
@@ -562,12 +498,12 @@ subroutine MPP_DO_GROUP_UPDATE_(group, domain, d_type)
                  if( 2*is-domain%x(1)%domain_data%begin.GT.domain%x(1)%domain_data%end+shift ) &
                       call mpp_error( FATAL, 'MPP_DO_UPDATE_V: folded-north BGRID_NE west edge ubound error.' )
                  do l=1,nvector
-                    ptr_fieldx = group%addrs_x(l)
-                    ptr_fieldy = group%addrs_y(l)
                     do k = 1,ksize
                        do i = domain%x(1)%domain_data%begin,is-1
-                          fieldx(i,j,k) = fieldx(2*is-i,j,k)
-                          fieldy(i,j,k) = fieldy(2*is-i,j,k)
+                          !fieldx(i,j,k) = fieldx(2*is-i,j,k)
+                          !fieldy(i,j,k) = fieldy(2*is-i,j,k)
+                          call SET_VALUE_ (group, FIELD_X, l, i, j, k, GET_VALUE_ (group, FIELD_X, l, 2*is-i, j, k))
+                          call SET_VALUE_ (group, FIELD_Y, l, i, j, k, GET_VALUE_ (group, FIELD_Y, l, 2*is-i, j, k))
                        end do
                     end do
                  end do
@@ -579,10 +515,10 @@ subroutine MPP_DO_GROUP_UPDATE_(group, domain, d_type)
                  if( 2*is-domain%x(1)%domain_data%begin-1.GT.domain%x(1)%domain_data%end ) &
                       call mpp_error( FATAL, 'MPP_DO_UPDATE_V: folded-north CGRID_NE west edge ubound error.' )
                  do l=1,nvector
-                    ptr_fieldy = group%addrs_y(l)
                     do k = 1,ksize
                        do i = isd,is-1
-                          fieldy(i,j,k) = fieldy(2*is-i-1,j,k)
+                          !fieldy(i,j,k) = fieldy(2*is-i-1,j,k)
+                          call SET_VALUE_ (group, FIELD_Y, l, i, j, k, GET_VALUE_ (group, FIELD_Y, l, 2*is-i-1, j, k))
                        end do
                     end do
                  end do
@@ -599,21 +535,21 @@ subroutine MPP_DO_GROUP_UPDATE_(group, domain, d_type)
               is = is + shift
               ie = ie + shift
               do l=1,nvector
-                 ptr_fieldx = group%addrs_x(l)
-                 ptr_fieldy = group%addrs_y(l)
                  do k = 1,ksize
                     do i = is,ie
-                       fieldx(i,j,k) = -fieldx(i,j,k)
-                       fieldy(i,j,k) = -fieldy(i,j,k)
+                       !fieldx(i,j,k) = -fieldx(i,j,k)
+                       !fieldy(i,j,k) = -fieldy(i,j,k)
+                       call SET_VALUE_ (group, FIELD_X, l, i, j, k, - GET_VALUE_ (group, FIELD_X, l, i, j, k))
+                       call SET_VALUE_ (group, FIELD_Y, l, i, j, k, - GET_VALUE_ (group, FIELD_Y, l, i, j, k))
                     end do
                  end do
               end do
            case(CGRID_NE)
               do l=1,nvector
-                 ptr_fieldy = group%addrs_y(l)
                  do k = 1,ksize
                     do i = is, ie
-                       fieldy(i,j,k) = -fieldy(i,j,k)
+                       !fieldy(i,j,k) = -fieldy(i,j,k)
+                       call SET_VALUE_ (group, FIELD_Y, l, i, j, k, - GET_VALUE_ (group, FIELD_Y, l, i, j, k))
                     end do
                  end do
               end do
@@ -648,15 +584,9 @@ subroutine MPP_START_GROUP_UPDATE_(group, domain, d_type, reuse_buffer)
   integer   :: n, l, m, i, j, k, buffer_start_pos, nk
   logical   :: reuse_buf_pos
   character(len=8)            :: text
+  type(c_ptr) :: stack_cptr !< Workaround for GFortran bug
 
-  MPP_TYPE_ :: buffer(size(mpp_domains_stack_nonblock(:)))
-  MPP_TYPE_ :: field (group%is_s:group%ie_s,group%js_s:group%je_s, group%ksize_s)
-  MPP_TYPE_ :: fieldx(group%is_x:group%ie_x,group%js_x:group%je_x, group%ksize_v)
-  MPP_TYPE_ :: fieldy(group%is_y:group%ie_y,group%js_y:group%je_y, group%ksize_v)
-  pointer( ptr, buffer )
-  pointer(ptr_field, field)
-  pointer(ptr_fieldx, fieldx)
-  pointer(ptr_fieldy, fieldy)
+  MPP_TYPE_, pointer :: buffer(:)
 
   nscalar = group%nscalar
   nvector = group%nvector
@@ -695,7 +625,8 @@ subroutine MPP_START_GROUP_UPDATE_(group, domain, d_type, reuse_buffer)
   nrecv = group%nrecv
   nsend = group%nsend
 
-  ptr = LOC(mpp_domains_stack_nonblock)
+  stack_cptr = c_loc(mpp_domains_stack_nonblock(1))
+  call c_f_pointer(stack_cptr, buffer, [size(mpp_domains_stack_nonblock)])
 
   ! Make sure it is not in the middle of the old version of non-blocking halo update.
   if(num_update>0) call mpp_error(FATAL, "MPP_START_GROUP_UPDATE: can not be called in the middle of "// &
@@ -751,14 +682,10 @@ subroutine MPP_COMPLETE_GROUP_UPDATE_(group, domain, d_type)
   integer   :: shift, gridtype, midpoint, flags_v
   integer   :: nunpack, rotation, buffer_start_pos, nk, isd
   logical   :: recv_y(8)
-  MPP_TYPE_ :: buffer(size(mpp_domains_stack_nonblock(:)))
-  MPP_TYPE_ :: field (group%is_s:group%ie_s,group%js_s:group%je_s, group%ksize_s)
-  MPP_TYPE_ :: fieldx(group%is_x:group%ie_x,group%js_x:group%je_x, group%ksize_v)
-  MPP_TYPE_ :: fieldy(group%is_y:group%ie_y,group%js_y:group%je_y, group%ksize_v)
-  pointer(ptr, buffer )
-  pointer(ptr_field, field)
-  pointer(ptr_fieldx, fieldx)
-  pointer(ptr_fieldy, fieldy)
+  type(c_ptr) :: stack_cptr !< Workaround for GFortran bug
+
+  MPP_TYPE_, pointer :: buffer(:)
+  MPP_TYPE_, parameter :: zero_ = 0.
 
   gridtype = group%gridtype
   flags_v = group%flags_v
@@ -772,7 +699,9 @@ subroutine MPP_COMPLETE_GROUP_UPDATE_(group, domain, d_type)
      ksize = group%ksize_v
   endif
   if(nvector > 0) recv_y = group%recv_y
-  ptr = LOC(mpp_domains_stack_nonblock)
+
+  stack_cptr = c_loc(mpp_domains_stack_nonblock(1))
+  call c_f_pointer(stack_cptr, buffer, [size(mpp_domains_stack_nonblock)])
 
   if(num_nonblock_group_update < 1) call mpp_error(FATAL, &
     'mpp_start_group_update must be called before calling mpp_end_group_update')
@@ -809,11 +738,9 @@ subroutine MPP_COMPLETE_GROUP_UPDATE_(group, domain, d_type)
            do i = is ,ie, midpoint
               if( domain%x(1)%domain_data%begin.LE.i .AND. i.LE. domain%x(1)%domain_data%end+shift )then
                  do l=1,nvector
-                    ptr_fieldx = group%addrs_x(l)
-                    ptr_fieldy = group%addrs_y(l)
                     do k = 1,ksize
-                       fieldx(i,j,k) = 0.
-                       fieldy(i,j,k) = 0.
+                       call SET_VALUE_ (group, FIELD_X, l, i, j, k, zero_)
+                       call SET_VALUE_ (group, FIELD_Y, l, i, j, k, zero_)
                     end do
                  end do
               end if
@@ -836,12 +763,12 @@ subroutine MPP_COMPLETE_GROUP_UPDATE_(group, domain, d_type)
                  if( 2*is-domain%x(1)%domain_data%begin.GT.domain%x(1)%domain_data%end+shift ) &
                       call mpp_error( FATAL, 'MPP_DO_UPDATE_V: folded-north BGRID_NE west edge ubound error.' )
                  do l=1,nvector
-                    ptr_fieldx = group%addrs_x(l)
-                    ptr_fieldy = group%addrs_y(l)
                     do k = 1,ksize
                        do i = domain%x(1)%domain_data%begin,is-1
-                          fieldx(i,j,k) = fieldx(2*is-i,j,k)
-                          fieldy(i,j,k) = fieldy(2*is-i,j,k)
+                          !fieldx(i,j,k) = fieldx(2*is-i,j,k)
+                          !fieldy(i,j,k) = fieldy(2*is-i,j,k)
+                          call SET_VALUE_ (group, FIELD_X, l, i, j, k, GET_VALUE_ (group, FIELD_X, l, 2*is-i, j, k))
+                          call SET_VALUE_ (group, FIELD_Y, l, i, j, k, GET_VALUE_ (group, FIELD_Y, l, 2*is-i, j, k))
                        end do
                     end do
                  end do
@@ -853,10 +780,10 @@ subroutine MPP_COMPLETE_GROUP_UPDATE_(group, domain, d_type)
                  if( 2*is-domain%x(1)%domain_data%begin-1.GT.domain%x(1)%domain_data%end ) &
                       call mpp_error( FATAL, 'MPP_DO_UPDATE_V: folded-north CGRID_NE west edge ubound error.' )
                  do l=1,nvector
-                    ptr_fieldy = group%addrs_y(l)
                     do k = 1,ksize
                        do i = isd,is-1
-                          fieldy(i,j,k) = fieldy(2*is-i-1,j,k)
+                          !fieldy(i,j,k) = fieldy(2*is-i-1,j,k)
+                          call SET_VALUE_ (group, FIELD_Y, l, i, j, k, GET_VALUE_ (group, FIELD_Y, l, 2*is-i-1, j, k))
                        end do
                     end do
                  end do
@@ -873,21 +800,21 @@ subroutine MPP_COMPLETE_GROUP_UPDATE_(group, domain, d_type)
               is = is + shift
               ie = ie + shift
               do l=1,nvector
-                 ptr_fieldx = group%addrs_x(l)
-                 ptr_fieldy = group%addrs_y(l)
                  do k = 1,ksize
                     do i = is,ie
-                       fieldx(i,j,k) = -fieldx(i,j,k)
-                       fieldy(i,j,k) = -fieldy(i,j,k)
+                       !fieldx(i,j,k) = -fieldx(i,j,k)
+                       !fieldy(i,j,k) = -fieldy(i,j,k)
+                       call SET_VALUE_ (group, FIELD_X, l, i, j, k, - GET_VALUE_ (group, FIELD_X, l, i, j, k))
+                       call SET_VALUE_ (group, FIELD_Y, l, i, j, k, - GET_VALUE_ (group, FIELD_Y, l, i, j, k))
                     end do
                  end do
               end do
            case(CGRID_NE)
               do l=1,nvector
-                 ptr_fieldy = group%addrs_y(l)
                  do k = 1,ksize
                     do i = is, ie
-                       fieldy(i,j,k) = -fieldy(i,j,k)
+                       !fieldy(i,j,k) = -fieldy(i,j,k)
+                       call SET_VALUE_ (group, FIELD_Y, l, i, j, k, - GET_VALUE_ (group, FIELD_Y, l, i, j, k))
                     end do
                  end do
               end do
@@ -911,114 +838,276 @@ subroutine MPP_COMPLETE_GROUP_UPDATE_(group, domain, d_type)
 
 end subroutine MPP_COMPLETE_GROUP_UPDATE_
 
-subroutine MPP_RESET_GROUP_UPDATE_FIELD_2D_(group, field)
-  type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,                   intent(in)    :: field(:,:)
-
-  group%reset_index_s = group%reset_index_s + 1
-
-  if(group%reset_index_s > group%nscalar) &
-     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_2D_: group%reset_index_s > group%nscalar")
-  if(size(field,1) .NE. group%isize_s .OR. size(field,2) .NE. group%jsize_s .OR. group%ksize_s .NE. 1) &
-     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_2D_: size of field does not match the size stored in group")
-
-  group%addrs_s(group%reset_index_s) = LOC(field)
-
-end subroutine MPP_RESET_GROUP_UPDATE_FIELD_2D_
-
-subroutine MPP_RESET_GROUP_UPDATE_FIELD_3D_(group, field)
-  type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,                   intent(in)    :: field(:,:,:)
-
-  group%reset_index_s = group%reset_index_s + 1
-
-  if(group%reset_index_s > group%nscalar) &
-     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_3D_: group%reset_index_s > group%nscalar")
-  if(size(field,1) .NE. group%isize_s .OR. size(field,2) .NE. group%jsize_s .OR. size(field,3) .NE. group%ksize_s) &
-     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_3D_: size of field does not match the size stored in group")
-
-  group%addrs_s(group%reset_index_s) = LOC(field)
-
-end subroutine MPP_RESET_GROUP_UPDATE_FIELD_3D_
-
-subroutine MPP_RESET_GROUP_UPDATE_FIELD_4D_(group, field)
+subroutine MPP_RESET_GROUP_UPDATE_FIELD_(group, field)
   type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,                   intent(in)    :: field(:,:,:,:)
+  MPP_TYPE_, target,           intent(in)    :: field(..)
+  integer :: isize, jsize, ksize
 
   group%reset_index_s = group%reset_index_s + 1
 
   if(group%reset_index_s > group%nscalar) &
-     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_4D_: group%reset_index_s > group%nscalar")
-  if(size(field,1) .NE. group%isize_s .OR. size(field,2) .NE. group%jsize_s .OR. &
-              size(field,3)*size(field,4) .NE. group%ksize_s) &
-     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_4D_: size of field does not match the size stored in group")
+     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_: group%reset_index_s > group%nscalar")
 
-  group%addrs_s(group%reset_index_s) = LOC(field)
+  isize = size(field, group%ix_s)
+  jsize = size(field, group%iy_s)
+  ksize = size(field) / (isize*jsize)
 
-end subroutine MPP_RESET_GROUP_UPDATE_FIELD_4D_
+  if(isize.ne.group%isize_s .or. jsize.ne.group%jsize_s .or. ksize.ne.group%ksize_s) then
+    call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_: size of field does not match the size stored in group")
+  endif
 
+  group%addrs_s(group%reset_index_s) = c_loc(field)
+end subroutine MPP_RESET_GROUP_UPDATE_FIELD_
 
-subroutine MPP_RESET_GROUP_UPDATE_FIELD_2D_V_(group, fieldx, fieldy)
+subroutine MPP_RESET_GROUP_UPDATE_FIELD_V_(group, fieldx, fieldy)
   type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,                   intent(in)    :: fieldx(:,:), fieldy(:,:)
+  MPP_TYPE_, target,           intent(in)    :: fieldx(..), fieldy(..)
+  integer :: isize_x, jsize_x, ksize_x
+  integer :: isize_y, jsize_y, ksize_y
 
   group%reset_index_v = group%reset_index_v + 1
 
   if(group%reset_index_v > group%nvector) &
-     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_2D_V_: group%reset_index_v > group%nvector")
-  if(size(fieldx,1) .NE. group%isize_x .OR. size(fieldx,2) .NE. group%jsize_x .OR. group%ksize_v .NE. 1) &
-     call mpp_error(FATAL, &
-                    &  "MPP_RESET_GROUP_UPDATE_FIELD_2D_V_: size of fieldx does not match the size stored in group")
-  if(size(fieldy,1) .NE. group%isize_y .OR. size(fieldy,2) .NE. group%jsize_y ) &
-     call mpp_error(FATAL, &
-                    &  "MPP_RESET_GROUP_UPDATE_FIELD_2D_V_: size of fieldy does not match the size stored in group")
-
-  group%addrs_x(group%reset_index_v) = LOC(fieldx)
-  group%addrs_y(group%reset_index_v) = LOC(fieldy)
+     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_V_: group%reset_index_v > group%nvector")
 
-end subroutine MPP_RESET_GROUP_UPDATE_FIELD_2D_V_
+  isize_x = size(fieldx, group%ix_v)
+  jsize_x = size(fieldx, group%iy_v)
+  ksize_x = size(fieldx) / (isize_x*jsize_x)
 
+  isize_y = size(fieldy, group%ix_v)
+  jsize_y = size(fieldy, group%iy_v)
+  ksize_y = size(fieldy) / (isize_y*jsize_y)
 
-subroutine MPP_RESET_GROUP_UPDATE_FIELD_3D_V_(group, fieldx, fieldy)
-  type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,                   intent(in)    :: fieldx(:,:,:), fieldy(:,:,:)
-
-  group%reset_index_v = group%reset_index_v + 1
-
-  if(group%reset_index_v > group%nvector) &
-     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_3D_V_: group%reset_index_v > group%nvector")
-  if(size(fieldx,1) .NE. group%isize_x .OR. size(fieldx,2) .NE. group%jsize_x .OR. size(fieldx,3) .NE. group%ksize_v) &
-     call mpp_error(FATAL, &
-                    &  "MPP_RESET_GROUP_UPDATE_FIELD_3D_V_: size of fieldx does not match the size stored in group")
-  if(size(fieldy,1) .NE. group%isize_y .OR. size(fieldy,2) .NE. group%jsize_y .OR. size(fieldy,3) .NE. group%ksize_v) &
-     call mpp_error(FATAL, &
-                    &  "MPP_RESET_GROUP_UPDATE_FIELD_3D_V_: size of fieldy does not match the size stored in group")
+  if(isize_x.ne.group%isize_x .or. jsize_x .ne. group%jsize_x .or. ksize_x .ne. group%ksize_v) then
+    call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_V_: size of fieldx does not match the size stored in group")
+  endif
 
-  group%addrs_x(group%reset_index_v) = LOC(fieldx)
-  group%addrs_y(group%reset_index_v) = LOC(fieldy)
+  if(isize_y.ne.group%isize_y .or. jsize_y .ne. group%jsize_y .or. ksize_y .ne. group%ksize_v) then
+    call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_V_: size of fieldy does not match the size stored in group")
+  endif
 
-end subroutine MPP_RESET_GROUP_UPDATE_FIELD_3D_V_
+  group%addrs_x(group%reset_index_v) = c_loc(fieldx)
+  group%addrs_y(group%reset_index_v) = c_loc(fieldy)
+end subroutine MPP_RESET_GROUP_UPDATE_FIELD_V_
+
+function GET_VALUE_ (group, field_type, l, i, j, k) result(res)
+  type(mpp_group_update_type), intent(in), target :: group
+  integer, intent(in) :: field_type !< FIELD_S, FIELD_X, or FIELD_Y
+  integer, intent(in) :: l !< Index of the field within group
+  integer, intent(in) :: i, j, k !< i,j,k indices of the value to retrieve
+  MPP_TYPE_ :: res
+
+  integer, allocatable, dimension(:) :: lb, n, indx
+  MPP_TYPE_, pointer :: field2d(:,:)
+  MPP_TYPE_, pointer :: field3d(:,:,:)
+  MPP_TYPE_, pointer :: field4d(:,:,:,:)
+
+  integer :: ix, iy !< Indices of the domain-decomposed dimensions
+  integer :: is, ie, js, je !< Starting and ending indices of the x and y dimensions
+  integer :: nd, m, kp, d
+
+  select case (field_type)
+    case (FIELD_S)
+      ix = group%ix_s
+      iy = group%iy_s
+      is = group%is_s
+      ie = group%ie_s
+      js = group%js_s
+      je = group%je_s
+
+      nd = size(group%shape_s)
+      allocate (lb(nd), n(nd), indx(nd))
+
+      lb = 1
+      n = group%shape_s
+
+      select case(nd)
+        case (2)
+          call c_f_pointer(group%addrs_s(l), field2d, group%shape_s)
+        case (3)
+          call c_f_pointer(group%addrs_s(l), field3d, group%shape_s)
+        case (4)
+          call c_f_pointer(group%addrs_s(l), field4d, group%shape_s)
+      end select
+    case (FIELD_X)
+      ix = group%ix_v
+      iy = group%iy_v
+      is = group%is_x
+      ie = group%ie_x
+      js = group%js_x
+      je = group%je_x
+
+      nd = size(group%shape_x)
+      allocate (lb(nd), n(nd), indx(nd))
+
+      lb = 1
+      n = group%shape_x
+
+      select case(nd)
+        case (2)
+          call c_f_pointer(group%addrs_x(l), field2d, group%shape_x)
+        case (3)
+          call c_f_pointer(group%addrs_x(l), field3d, group%shape_x)
+        case (4)
+          call c_f_pointer(group%addrs_x(l), field4d, group%shape_x)
+      end select
+    case (FIELD_Y)
+      ix = group%ix_v
+      iy = group%iy_v
+      is = group%is_y
+      ie = group%ie_y
+      js = group%js_y
+      je = group%je_y
+
+      nd = size(group%shape_y)
+      allocate (lb(nd), n(nd), indx(nd))
+
+      lb = 1
+      n = group%shape_y
+
+      select case(nd)
+        case (2)
+          call c_f_pointer(group%addrs_y(l), field2d, group%shape_y)
+        case (3)
+          call c_f_pointer(group%addrs_y(l), field3d, group%shape_y)
+        case (4)
+          call c_f_pointer(group%addrs_y(l), field4d, group%shape_y)
+      end select
+  end select
+
+  lb(ix) = i - is + 1
+  lb(iy) = j - js + 1
+
+  n(ix) = 1
+  n(iy) = 1
+
+  kp = k - 1
+  do d=nd,1,-1
+    m = product(n(1:d-1)) ! Product of all lower dimension sizes
+    indx(d) = kp / m
+    kp = kp - indx(d)*m
+  enddo
 
+  indx = lb + indx
+
+  select case (nd)
+    case (2)
+      res = field2d(indx(1), indx(2))
+    case (3)
+      res = field3d(indx(1), indx(2), indx(3))
+    case (4)
+      res = field4d(indx(1), indx(2), indx(3), indx(4))
+  end select
+end function GET_VALUE_
+
+subroutine SET_VALUE_ (group, field_type, l, i, j, k, x)
+  type(mpp_group_update_type), intent(in), target :: group
+  integer, intent(in) :: field_type !< FIELD_S, FIELD_X, or FIELD_Y
+  integer, intent(in) :: l !< Index of the field within group
+  integer, intent(in) :: i, j, k !< i,j,k indices of the value to retrieve
+  MPP_TYPE_, intent(in) :: x
+
+  integer, allocatable, dimension(:) :: lb, n, indx
+  MPP_TYPE_, pointer :: field2d(:,:)
+  MPP_TYPE_, pointer :: field3d(:,:,:)
+  MPP_TYPE_, pointer :: field4d(:,:,:,:)
+
+  integer :: ix, iy !< Indices of the domain-decomposed dimensions
+  integer :: is, ie, js, je !< Starting and ending indices of the x and y dimensions
+  integer :: nd, m, kp, d
+
+  select case (field_type)
+    case (FIELD_S)
+      ix = group%ix_s
+      iy = group%iy_s
+      is = group%is_s
+      ie = group%ie_s
+      js = group%js_s
+      je = group%je_s
+
+      nd = size(group%shape_s)
+      allocate (lb(nd), n(nd), indx(nd))
+
+      lb = 1
+      n = group%shape_s
+
+      select case(nd)
+        case (2)
+          call c_f_pointer(group%addrs_s(l), field2d, group%shape_s)
+        case (3)
+          call c_f_pointer(group%addrs_s(l), field3d, group%shape_s)
+        case (4)
+          call c_f_pointer(group%addrs_s(l), field4d, group%shape_s)
+      end select
+    case (FIELD_X)
+      ix = group%ix_v
+      iy = group%iy_v
+      is = group%is_x
+      ie = group%ie_x
+      js = group%js_x
+      je = group%je_x
+
+      nd = size(group%shape_x)
+      allocate (lb(nd), n(nd), indx(nd))
+
+      lb = 1
+      n = group%shape_x
+
+      select case(nd)
+        case (2)
+          call c_f_pointer(group%addrs_x(l), field2d, group%shape_x)
+        case (3)
+          call c_f_pointer(group%addrs_x(l), field3d, group%shape_x)
+        case (4)
+          call c_f_pointer(group%addrs_x(l), field4d, group%shape_x)
+      end select
+    case (FIELD_Y)
+      ix = group%ix_v
+      iy = group%iy_v
+      is = group%is_y
+      ie = group%ie_y
+      js = group%js_y
+      je = group%je_y
+
+      nd = size(group%shape_y)
+      allocate (lb(nd), n(nd), indx(nd))
+
+      lb = 1
+      n = group%shape_y
+
+      select case(nd)
+        case (2)
+          call c_f_pointer(group%addrs_y(l), field2d, group%shape_y)
+        case (3)
+          call c_f_pointer(group%addrs_y(l), field3d, group%shape_y)
+        case (4)
+          call c_f_pointer(group%addrs_y(l), field4d, group%shape_y)
+      end select
+  end select
+
+  lb(ix) = i - is + 1
+  lb(iy) = j - js + 1
+
+  n(ix) = 1
+  n(iy) = 1
+
+  kp = k - 1
+  do d=nd,1,-1
+    m = product(n(1:d-1)) ! Product of all lower dimension sizes
+    indx(d) = kp / m
+    kp = kp - indx(d)*m
+  enddo
 
-subroutine MPP_RESET_GROUP_UPDATE_FIELD_4D_V_(group, fieldx, fieldy)
-  type(mpp_group_update_type), intent(inout) :: group
-  MPP_TYPE_,                   intent(in)    :: fieldx(:,:,:,:), fieldy(:,:,:,:)
+  indx = lb + indx
 
-  group%reset_index_v = group%reset_index_v + 1
+  select case (nd)
+    case (2)
+      field2d(indx(1), indx(2)) = x
+    case (3)
+      field3d(indx(1), indx(2), indx(3)) = x
+    case (4)
+      field4d(indx(1), indx(2), indx(3), indx(4)) = x
+  end select
+end subroutine SET_VALUE_
 
-  if(group%reset_index_v > group%nvector) &
-     call mpp_error(FATAL, "MPP_RESET_GROUP_UPDATE_FIELD_4D_V_: group%reset_index_v > group%nvector")
-  if(size(fieldx,1) .NE. group%isize_x .OR. size(fieldx,2) .NE. group%jsize_x .OR. &
-              size(fieldx,3)*size(fieldx,4) .NE. group%ksize_v) &
-     call mpp_error(FATAL, &
-                    &  "MPP_RESET_GROUP_UPDATE_FIELD_4D_V_: size of fieldx does not match the size stored in group")
-  if(size(fieldy,1) .NE. group%isize_y .OR. size(fieldy,2) .NE. group%jsize_y .OR. &
-              size(fieldy,3)*size(fieldy,4) .NE. group%ksize_v) &
-     call mpp_error(FATAL, &
-                    &  "MPP_RESET_GROUP_UPDATE_FIELD_4D_V_: size of fieldy does not match the size stored in group")
-
-  group%addrs_x(group%reset_index_v) = LOC(fieldx)
-  group%addrs_y(group%reset_index_v) = LOC(fieldy)
-
-end subroutine MPP_RESET_GROUP_UPDATE_FIELD_4D_V_
 !> @}
diff --git a/mpp/mpp_domains.F90 b/mpp/mpp_domains.F90
index 722430f1b4..2990e084d6 100644
--- a/mpp/mpp_domains.F90
+++ b/mpp/mpp_domains.F90
@@ -96,6 +96,7 @@ module mpp_domains_mod
   use gfdl_nompi_f08
 #endif
 
+  use iso_c_binding,          only : c_f_pointer, c_loc, c_ptr
   use mpp_parameter_mod,      only : MPP_DEBUG, MPP_VERBOSE, MPP_DOMAIN_TIME
   use mpp_parameter_mod,      only : GLOBAL_DATA_DOMAIN, CYCLIC_GLOBAL_DOMAIN, GLOBAL,CYCLIC
   use mpp_parameter_mod,      only : AGRID, BGRID_SW, BGRID_NE, CGRID_NE, CGRID_SW, DGRID_NE, DGRID_SW
@@ -592,6 +593,11 @@ module mpp_domains_mod
      integer            :: is_s=0, ie_s=0, js_s=0, je_s=0
      integer            :: is_x=0, ie_x=0, js_x=0, je_x=0
      integer            :: is_y=0, ie_y=0, js_y=0, je_y=0
+     integer            :: ix_s=1, iy_s=2 !< Domain-decomposed dimensions of scalar arrays
+     integer            :: ix_v=1, iy_v=2 !< Domain-decomposed dimensions of vector arrays
+     integer, allocatable, dimension(:) :: shape_s !< Shapes of scalar arrays
+     integer, allocatable, dimension(:) :: shape_x !< Shapes of vector arrays (x-component)
+     integer, allocatable, dimension(:) :: shape_y !< Shapes of vector arrays (y-component)
      integer            :: nrecv=0, nsend=0
      integer            :: npack=0, nunpack=0
      integer            :: reset_index_s = 0
@@ -619,9 +625,9 @@ module mpp_domains_mod
      integer            :: unpack_ie(MAXOVERLAP)
      integer            :: unpack_js(MAXOVERLAP)
      integer            :: unpack_je(MAXOVERLAP)
-     integer(i8_kind)   :: addrs_s(MAX_DOMAIN_FIELDS)
-     integer(i8_kind)   :: addrs_x(MAX_DOMAIN_FIELDS)
-     integer(i8_kind)   :: addrs_y(MAX_DOMAIN_FIELDS)
+     type(c_ptr)        :: addrs_s(MAX_DOMAIN_FIELDS)
+     type(c_ptr)        :: addrs_x(MAX_DOMAIN_FIELDS)
+     type(c_ptr)        :: addrs_y(MAX_DOMAIN_FIELDS)
      integer            :: buffer_start_pos = -1
      type(mpi_request)  :: request_send(MAX_REQUEST)
      type(mpi_request)  :: request_recv(MAX_REQUEST)
@@ -1314,18 +1320,10 @@ module mpp_domains_mod
   !> @param
   !> @ingroup mpp_domains_mod
   interface mpp_create_group_update
-     module procedure mpp_create_group_update_r4_2d
-     module procedure mpp_create_group_update_r4_3d
-     module procedure mpp_create_group_update_r4_4d
-     module procedure mpp_create_group_update_r4_2dv
-     module procedure mpp_create_group_update_r4_3dv
-     module procedure mpp_create_group_update_r4_4dv
-     module procedure mpp_create_group_update_r8_2d
-     module procedure mpp_create_group_update_r8_3d
-     module procedure mpp_create_group_update_r8_4d
-     module procedure mpp_create_group_update_r8_2dv
-     module procedure mpp_create_group_update_r8_3dv
-     module procedure mpp_create_group_update_r8_4dv
+     module procedure mpp_create_group_update_r4
+     module procedure mpp_create_group_update_r4_v
+     module procedure mpp_create_group_update_r8
+     module procedure mpp_create_group_update_r8_v
   end interface mpp_create_group_update
 
   !> @ingroup mpp_domains_mod
@@ -1360,18 +1358,10 @@ module mpp_domains_mod
 
   !> @ingroup mpp_domains_mod
   interface mpp_reset_group_update_field
-     module procedure mpp_reset_group_update_field_r4_2d
-     module procedure mpp_reset_group_update_field_r4_3d
-     module procedure mpp_reset_group_update_field_r4_4d
-     module procedure mpp_reset_group_update_field_r4_2dv
-     module procedure mpp_reset_group_update_field_r4_3dv
-     module procedure mpp_reset_group_update_field_r4_4dv
-     module procedure mpp_reset_group_update_field_r8_2d
-     module procedure mpp_reset_group_update_field_r8_3d
-     module procedure mpp_reset_group_update_field_r8_4d
-     module procedure mpp_reset_group_update_field_r8_2dv
-     module procedure mpp_reset_group_update_field_r8_3dv
-     module procedure mpp_reset_group_update_field_r8_4dv
+     module procedure mpp_reset_group_update_field_r4
+     module procedure mpp_reset_group_update_field_r4_v
+     module procedure mpp_reset_group_update_field_r8
+     module procedure mpp_reset_group_update_field_r8_v
   end interface mpp_reset_group_update_field
 
   !> Pass the data from coarse grid to fill the buffer to be ready to be interpolated
diff --git a/test_fms/mpp/include/group_update.inc b/test_fms/mpp/include/group_update.inc
index 575a0f028b..7061cacf98 100644
--- a/test_fms/mpp/include/group_update.inc
+++ b/test_fms/mpp/include/group_update.inc
@@ -60,6 +60,7 @@
     integer        :: i, j, k, l, n, shift
     integer        :: isc, iec, jsc, jec, isd, ied, jsd, jed
     integer        :: ism, iem, jsm, jem
+    integer        :: indx(3), ix, iy
 
     integer, allocatable, dimension(:) :: pe_start, pe_end, tile1, tile2
     integer, allocatable, dimension(:) :: istart1, iend1, jstart1, jend1
@@ -170,6 +171,11 @@
     call mpp_get_data_domain   ( domain, isd, ied, jsd, jed )
     call mpp_get_memory_domain   ( domain, ism, iem, jsm, jem )
 
+    indx(1:3) = [1, 2, 3]
+    call permute_arr(indx, p)
+    ix = findloc(indx, 1, dim=1)
+    iy = findloc(indx, 2, dim=1)
+
     mem%lb = [ism, jsm, 1]
     mem%ub = [iem, jem, nz]
     call mem%permute(p)
@@ -248,7 +254,7 @@
 
     !--- Test for partial direction update
     do l =1, num_fields
-       call mpp_create_group_update(group_update, a1(:,:,:,l), domain, flags=WUPDATE+SUPDATE)
+       call mpp_create_group_update(group_update, a1(:,:,:,l), domain, flags=WUPDATE+SUPDATE, xdim=ix, ydim=iy)
     end do
 
     do l = 1, num_fields
@@ -276,7 +282,8 @@
        y1 = zero
 
        do l =1, num_fields
-          call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=DGRID_NE)
+          call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=DGRID_NE, &
+                                       xdim=ix, ydim=iy)
        end do
 
        do l = 1, num_fields
@@ -308,8 +315,9 @@
     y1 = zero
 
     do l =1, num_fields
-       call mpp_create_group_update(group_update, a1(:,:,:,l), domain)
-       call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=CGRID_NE)
+       call mpp_create_group_update(group_update, a1(:,:,:,l), domain, xdim=ix, ydim=iy)
+       call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=CGRID_NE, &
+                                    xdim=ix, ydim=iy)
     end do
 
     do n = 1, num_iter
@@ -375,8 +383,9 @@
     !--- The following is to test overlapping start and complete
     if( num_fields > 1 ) then
        do l =1, num_fields
-          call mpp_create_group_update(update_list(l), a1(:,:,:,l), domain)
-          call mpp_create_group_update(update_list(l), x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=CGRID_NE)
+          call mpp_create_group_update(update_list(l), a1(:,:,:,l), domain, xdim=ix, ydim=iy)
+          call mpp_create_group_update(update_list(l), x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=CGRID_NE, &
+                                       xdim=ix, ydim=iy)
        end do
 
        do n = 1, num_iter
@@ -478,8 +487,9 @@
     allocate(y2(INDICES_MEM_IJSHIFT_, num_fields))
 
     do l =1, num_fields
-       call mpp_create_group_update(group_update, a1(:,:,:,l), domain, position=CORNER)
-       call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=BGRID_NE)
+       call mpp_create_group_update(group_update, a1(:,:,:,l), domain, position=CORNER, xdim=ix, ydim=iy)
+       call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=BGRID_NE, &
+                                    xdim=ix, ydim=iy)
     end do
 
     do n = 1, num_iter
@@ -569,7 +579,7 @@
     y2 = y1
 
     do l =1, num_fields
-       call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=AGRID)
+       call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=AGRID, xdim=ix, ydim=iy)
     end do
 
     do l = 1, num_fields
@@ -600,7 +610,8 @@
     y2 = y1
 
     do l =1, num_fields
-       call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=AGRID, flags=SCALAR_PAIR)
+       call mpp_create_group_update(group_update, x1(:,:,:,l), y1(:,:,:,l), domain, gridtype=AGRID, &
+                                    flags=SCALAR_PAIR, xdim=ix, ydim=iy)
     end do
 
     do l = 1, num_fields
diff --git a/test_fms/mpp/test_mpp_domains.F90 b/test_fms/mpp/test_mpp_domains.F90
index b9ba0f3b0a..e681173944 100644
--- a/test_fms/mpp/test_mpp_domains.F90
+++ b/test_fms/mpp/test_mpp_domains.F90
@@ -58,7 +58,7 @@ program test_mpp_domains
   use compare_data_checksums
   use test_domains_utility_mod
   use platform_mod
-  use fms_test_mod, only: permutable_indices_3d, factorial, arr_init, arr_compare
+  use fms_test_mod, only: permutable_indices_3d, permute_arr, factorial, arr_init, arr_compare
 
 
   implicit none