diff --git a/GeomIO/CMakeLists.txt b/GeomIO/CMakeLists.txt
index 10c45bc1de62..bdcab8003489 100644
--- a/GeomIO/CMakeLists.txt
+++ b/GeomIO/CMakeLists.txt
@@ -6,6 +6,7 @@ set(srcs
   Geom_PFIO.F90
   Grid_PFIO.F90
   GeomCatagorizer.F90
+  pFIOServerBounds.F90
   )
 
 esma_add_library(${this}
diff --git a/GeomIO/Grid_PFIO.F90 b/GeomIO/Grid_PFIO.F90
index 88933d46e2d2..c94975d79a82 100644
--- a/GeomIO/Grid_PFIO.F90
+++ b/GeomIO/Grid_PFIO.F90
@@ -3,9 +3,13 @@
 module mapl3g_GridPFIO
    use mapl_ErrorHandling
    use mapl3g_GeomPFIO
+   use mapl3g_SharedIO
    use ESMF
    use PFIO
    use MAPL_BaseMod
+   use MAPL_FieldPointerUtilities
+   use mapl3g_pFIOServerBounds
+   use, intrinsic :: iso_c_binding, only: c_ptr
    implicit none
    private
 
@@ -30,11 +34,14 @@ subroutine stage_data_to_file(this, bundle, filename, time_index, rc)
       character(len=ESMF_MAXSTR), allocatable :: field_names(:)
       type(ESMF_Field) :: field
       type(ArrayReference) :: ref
-      real, pointer :: ptr2d(:,:)
       integer, allocatable :: local_start(:), global_start(:), global_count(:)
+      type(c_ptr) :: address
+      integer :: type_kind
+      type(ESMF_TypeKind_Flag) :: tk
+      integer, allocatable :: element_count(:), new_element_count(:) 
 
       type(ESMF_Grid) :: grid
-      integer :: global_dim(3), i1, j1, in, jn
+      type(pFIOServerBounds) :: server_bounds
 
       collection_id = this%get_collection_id()
       call ESMF_FieldBundleGet(bundle, fieldCount=num_fields, _RC)
@@ -42,16 +49,21 @@ subroutine stage_data_to_file(this, bundle, filename, time_index, rc)
       call ESMF_FieldBundleGet(bundle, fieldNameList=field_names, _RC)
       do i=1,num_fields
          call ESMF_FieldBundleGet(bundle, field_names(i), field=field, _RC)
-         ! all this logic needs to be generalized
-         call ESMF_FieldGet(field, farrayPtr=ptr2d, _RC)
-         allocate(global_start, source=[1,1,time_index])
-         call ESMF_FieldGet(field, grid=grid, _RC)
-         call MAPL_GridGet(grid, globalCellCountPerDim=global_dim, _RC)
-         allocate(global_count, source=[global_dim(1),global_dim(2),1])
-         call MAPL_GridGetInterior(grid, i1, in, j1, jn)
-         allocate(local_start, source=[i1, j1,1])
-         ref = ArrayReference(ptr2d)
-         ! end generalization
+
+         element_count = FieldGetLocalElementCount(field, _RC)
+         call ESMF_FieldGet(field, grid=grid, typekind=tk,  _RC)
+
+         call server_bounds%initialize(grid, element_count, time_index=time_index, _RC)
+         global_start = server_bounds%get_global_start()
+         global_count = server_bounds%get_global_count()
+         local_start = server_bounds%get_local_start()
+         
+         ! generate array reference
+         call FieldGetCptr(field, address, _RC)
+         type_kind = esmf_to_pfio_type(tk, _RC)
+         new_element_count = server_bounds%get_file_shape()
+         ref = ArrayReference(address, type_kind, new_element_count)
+
          call o_clients%collective_stage_data(collection_id,filename, trim(field_names(i)), &
               ref, start=local_start, global_start=global_start, global_count=global_count)
       enddo
diff --git a/GeomIO/SharedIO.F90 b/GeomIO/SharedIO.F90
index 250f2c7833e1..7b0e3fe4b44e 100644
--- a/GeomIO/SharedIO.F90
+++ b/GeomIO/SharedIO.F90
@@ -5,6 +5,7 @@ module mapl3g_SharedIO
    use pfio
    use gFTL2_StringVector
    use mapl3g_geom_mgr
+   use MAPL_BaseMod
 
    implicit none
 
@@ -13,6 +14,7 @@ module mapl3g_SharedIO
    public get_mapl_geom
    public create_time_variable
    public bundle_to_metadata
+   public esmf_to_pfio_type
 
    contains
 
@@ -86,7 +88,6 @@ subroutine add_variable(metadata, field,  rc)
       mapl_geom => get_mapl_geom(esmfgeom, _RC)
       grid_variables = mapl_geom%get_gridded_dims()
       dims = string_vec_to_comma_sep(grid_variables)
-      dims = 'lon,lat'
       call ESMF_FieldGet(field, name=fname, typekind = typekind, _RC)
       ! add vertical dimension
 
diff --git a/GeomIO/pFIOServerBounds.F90 b/GeomIO/pFIOServerBounds.F90
new file mode 100644
index 000000000000..b8fad0db644a
--- /dev/null
+++ b/GeomIO/pFIOServerBounds.F90
@@ -0,0 +1,118 @@
+#include "MAPL_Generic.h"
+module mapl3g_pFIOServerBounds
+   use mapl_ErrorHandlingMod
+   use esmf
+   use pfio
+   use gFTL2_StringVector
+   use MAPL_BaseMod
+
+   implicit none
+   private
+
+   public :: pFIOServerBounds
+
+   integer, parameter :: grid_dims = 2
+
+   type :: pFIOServerBounds
+      private
+      integer, allocatable :: local_start(:)
+      integer, allocatable :: global_start(:)
+      integer, allocatable :: global_count(:)
+      integer, allocatable :: file_shape(:)
+   contains
+      procedure :: initialize
+      procedure :: get_local_start
+      procedure :: get_global_start
+      procedure :: get_global_count
+      procedure :: get_file_shape
+   end type pFIOServerBounds
+
+   contains
+
+   function get_local_start(this) result(local_start)
+      integer, allocatable :: local_start(:)
+      class(pFIOServerBounds), intent(in) :: this
+      local_start =this%local_start
+   end function get_local_start
+
+   function get_global_start(this) result(global_start)
+      integer, allocatable :: global_start(:)
+      class(pFIOServerBounds), intent(in) :: this
+      global_start =this%global_start
+   end function get_global_start
+
+   function get_global_count(this) result(global_count)
+      integer, allocatable :: global_count(:)
+      class(pFIOServerBounds), intent(in) :: this
+      global_count =this%global_count
+   end function get_global_count
+
+   function get_file_shape(this) result(file_shape)
+      integer, allocatable :: file_shape(:)
+      class(pFIOServerBounds), intent(in) :: this
+      file_shape =this%file_shape
+   end function get_file_shape
+
+   subroutine initialize(this, grid, field_shape, time_index, rc)
+      class(pFIOServerBounds), intent(inout) :: this
+      type(ESMF_Grid), intent(in) :: grid
+      integer, intent(in) :: field_shape(:)
+      integer, intent(in), optional :: time_index
+      integer, intent(out), optional :: rc
+
+      integer :: status, tile_count, n_dims, tm, global_dim(3)
+      integer :: i1, in, j1, jn, tile, extra_file_dim, file_dims, new_grid_dims
+
+      call ESMF_GridGet(grid, tileCount=tile_count, _RC)
+      call MAPL_GridGetInterior(grid, i1,in, j1, jn)
+      call MAPL_GridGet(grid, globalCellCountPerDim=global_dim, _RC)
+      n_dims = size(field_shape)
+
+      tm = 0
+      if (present(time_index))  tm = 1
+
+      extra_file_dim = 0
+      if (tile_count == 6) extra_file_dim = 1
+
+      new_grid_dims = grid_dims + extra_file_dim
+      file_dims = n_dims + extra_file_dim
+
+      allocate(this%file_shape(file_dims))
+      allocate(this%global_start(file_dims+tm))
+      allocate(this%global_count(file_dims+tm))
+      allocate(this%local_start(file_dims+tm))
+
+      this%file_shape(new_grid_dims+1:file_dims) = field_shape(grid_dims+1:n_dims)
+
+      this%global_start(1:file_dims) = 1
+      if(present(time_index)) this%global_start(file_dims+1) = time_index
+
+      this%global_count(new_grid_dims+1:file_dims) = field_shape(grid_dims+1:n_dims)
+      if (present(time_index)) this%global_count(file_dims+1) = 1 
+
+      this%local_start = 1
+
+      select case (tile_count)
+      case (6) ! Assume cubed-sphere
+
+         tile = 1 + (j1-1)/global_dim(1)
+         this%file_shape(1:new_grid_dims) = [field_shape(1), field_shape(2) ,1]
+         this%global_count(1:new_grid_dims) =[global_dim(1), global_dim(1), tile_count]
+         this%local_start(1:new_grid_dims) = [i1, j1-(tile-1)*global_dim(1), tile]
+
+      case (1) 
+         
+         this%file_shape(1:new_grid_dims) = [field_shape(1), field_shape(2)]
+         this%global_count(1:new_grid_dims) = [global_dim(1), global_dim(2)]
+         this%local_start(1:new_grid_dims) = [i1,j1]
+
+      case default
+         _FAIL("unsupported grid")
+      end select
+
+      _RETURN(_SUCCESS)
+
+   end subroutine initialize
+      
+end module mapl3g_pFIOServerBounds
+
diff --git a/pfio/ArrayReference.F90 b/pfio/ArrayReference.F90
index 67a9635ea132..92b149608957 100644
--- a/pfio/ArrayReference.F90
+++ b/pfio/ArrayReference.F90
@@ -3,7 +3,7 @@
 
 module pFIO_ArrayReferenceMod
    use, intrinsic :: iso_c_binding, only: C_NULL_PTR
-   use, intrinsic :: iso_c_binding, only: c_loc
+   use, intrinsic :: iso_c_binding, only: c_loc, c_ptr
    use, intrinsic :: iso_fortran_env, only: INT32
    use, intrinsic :: iso_fortran_env, only: INT64
    use, intrinsic :: iso_fortran_env, only: REAL32
@@ -25,16 +25,26 @@ module pFIO_ArrayReferenceMod
    end type ArrayReference
 
    interface ArrayReference
-      module procedure new_ArrayReference_0d
-      module procedure new_ArrayReference_1d
-      module procedure new_ArrayReference_2d
-      module procedure new_ArrayReference_3d
-      module procedure new_ArrayReference_4d
-      module procedure new_ArrayReference_5d
+      procedure new_ArrayReference_from_param
+      procedure new_ArrayReference_0d
+      procedure new_ArrayReference_1d
+      procedure new_ArrayReference_2d
+      procedure new_ArrayReference_3d
+      procedure new_ArrayReference_4d
+      procedure new_ArrayReference_5d
    end interface ArrayReference
 
 contains
 
+   function new_ArrayReference_from_param(in_c_loc, in_kind, in_shape) result(reference)
+      type (ArrayReference) :: reference
+      type(c_ptr), intent(in) :: in_c_loc
+      integer, intent(in) :: in_kind
+      integer, intent(in) :: in_shape(:)
+      reference%base_address = in_c_loc
+      reference%shape = in_shape
+      reference%type_kind = in_kind
+   end function
 
    function new_ArrayReference_0d(scalar, rc) result(reference)
       type (ArrayReference) :: reference