mpp_domains_util.inc

! -*-f90-*-
 
!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************
 
!> @file
!> @brief Utility routines for getting and setting values in @ref mpp_domains_mod
 
!> @ingroup mpp_domains_mod
!> @{
 
  !> @brief Set user stack size.
  !!
  !> This sets the size of an array that is used for internal storage by
  !! <TT>mpp_domains</TT>. This array is used, for instance, to buffer the
  !! data sent and received in halo updates.<br>
  !! This call has implied global synchronization. It should be
  !! placed somewhere where all PEs can call it.
  subroutine mpp_domains_set_stack_size(n)
    !set the mpp_domains_stack variable to be at least n LONG words long
    integer, intent(in) :: n
    character(len=8) :: text
 
    if( n.LE.mpp_domains_stack_size )return
    if( allocated(mpp_domains_stack) )deallocate(mpp_domains_stack)
    allocate( mpp_domains_stack(n) )
    if( allocated(mpp_domains_stack_nonblock) )deallocate(mpp_domains_stack_nonblock)
    allocate( mpp_domains_stack_nonblock(n) )
 
    mpp_domains_stack_size = n
    write( text,'(i8)' )n
    if( mpp_pe().EQ.mpp_root_pe() )call mpp_error( note, 'MPP_DOMAINS_SET_STACK_SIZE: stack size set to '//text//'.' )
 
    return
  end subroutine mpp_domains_set_stack_size
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                             !
  !                MPP_DOMAINS: overloaded operators (==, /=)                   !
  !                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  function mpp_domain1d_eq( a, b )
    logical                    :: mpp_domain1d_eq
    type(domain1d), intent(in) :: a, b
 
    mpp_domain1d_eq = ( a%compute%begin.EQ.b%compute%begin .AND. &
         a%compute%end  .EQ.b%compute%end   .AND. &
         a%domain_data%begin   .EQ.b%domain_data%begin    .AND. &
         a%domain_data%end     .EQ.b%domain_data%end      .AND. &
         a%global%begin .EQ.b%global%begin  .AND. &
         a%global%end   .EQ.b%global%end    )
    !compare pelists
    !      if( mpp_domain1D_eq )mpp_domain1D_eq = ASSOCIATED(a%list) .AND. ASSOCIATED(b%list)
    !      if( mpp_domain1D_eq )mpp_domain1D_eq = size(a%list(:)).EQ.size(b%list(:))
    !      if( mpp_domain1D_eq )mpp_domain1D_eq = ALL(a%list%pe.EQ.b%list%pe)
 
    return
  end function mpp_domain1d_eq
 
  function mpp_domain1d_ne( a, b )
    logical                    :: mpp_domain1d_ne
    type(domain1d), intent(in) :: a, b
 
    mpp_domain1d_ne = .NOT. ( a.EQ.b )
    return
  end function mpp_domain1d_ne
 
  function mpp_domain2d_eq( a, b )
    logical                    :: mpp_domain2d_eq
    type(domain2d), intent(in) :: a, b
    integer                    :: nt, n
 
    mpp_domain2d_eq = size(a%x(:)) .EQ. size(b%x(:))
    nt = size(a%x(:))
    do n = 1, nt
       if(mpp_domain2d_eq) mpp_domain2d_eq = a%x(n).EQ.b%x(n) .AND. a%y(n).EQ.b%y(n)
    end do
 
    if( mpp_domain2d_eq .AND. ((a%pe.EQ.null_pe).OR.(b%pe.EQ.null_pe)) )return !NULL_DOMAIN2D
    !compare pelists
    if( mpp_domain2d_eq )mpp_domain2d_eq = ASSOCIATED(a%list) .AND. ASSOCIATED(b%list)
    if( mpp_domain2d_eq )mpp_domain2d_eq = size(a%list(:)).EQ.size(b%list(:))
    if( mpp_domain2d_eq )mpp_domain2d_eq = all(a%list%pe.EQ.b%list%pe)
    if( mpp_domain2d_eq )mpp_domain2d_eq = all(a%io_layout .EQ. b%io_layout)
    if( mpp_domain2d_eq )mpp_domain2d_eq = a%symmetry .eqv. b%symmetry
 
    return
  end function mpp_domain2d_eq
 
  !#####################################################################
 
  function mpp_domain2d_ne( a, b )
    logical                    :: mpp_domain2d_ne
    type(domain2d), intent(in) :: a, b
 
    mpp_domain2d_ne = .NOT. ( a.EQ.b )
    return
  end function mpp_domain2d_ne
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                             !
  !     MPP_GET and SET routiness: retrieve various components of domains       !
  !                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  subroutine mpp_get_compute_domain1d( domain, begin, end, size, max_size, is_global )
    type(domain1d),     intent(in) :: domain
    integer, intent(out), optional :: begin, end, size, max_size
    logical, intent(out), optional :: is_global
 
    if( PRESENT(begin)     )begin     = domain%compute%begin
    if( PRESENT(end)       )end       = domain%compute%end
    if( PRESENT(size)      )size      = domain%compute%size
    if( PRESENT(max_size)  )max_size  = domain%compute%max_size
    if( PRESENT(is_global) )is_global = domain%compute%is_global
    return
  end subroutine mpp_get_compute_domain1d
 
  !#####################################################################
  subroutine mpp_get_data_domain1d( domain, begin, end, size, max_size, is_global )
    type(domain1d),     intent(in) :: domain
    integer, intent(out), optional :: begin, end, size, max_size
    logical, intent(out), optional :: is_global
 
    if( PRESENT(begin)     )begin     = domain%domain_data%begin
    if( PRESENT(end)       )end       = domain%domain_data%end
    if( PRESENT(size)      )size      = domain%domain_data%size
    if( PRESENT(max_size)  )max_size  = domain%domain_data%max_size
    if( PRESENT(is_global) )is_global = domain%domain_data%is_global
    return
  end subroutine mpp_get_data_domain1d
 
  !#####################################################################
  subroutine mpp_get_global_domain1d( domain, begin, end, size, max_size )
    type(domain1d),     intent(in) :: domain
    integer, intent(out), optional :: begin, end, size, max_size
 
    if( PRESENT(begin)    )begin    = domain%global%begin
    if( PRESENT(end)      )end      = domain%global%end
    if( PRESENT(size)     )size     = domain%global%size
    if( PRESENT(max_size) )max_size = domain%global%max_size
    return
  end subroutine mpp_get_global_domain1d
 
  !#####################################################################
  subroutine mpp_get_memory_domain1d( domain, begin, end, size, max_size, is_global )
    type(domain1d),     intent(in) :: domain
    integer, intent(out), optional :: begin, end, size, max_size
    logical, intent(out), optional :: is_global
 
    if( PRESENT(begin)     )begin     = domain%memory%begin
    if( PRESENT(end)       )end       = domain%memory%end
    if( PRESENT(size)      )size      = domain%memory%size
    if( PRESENT(max_size)  )max_size  = domain%memory%max_size
    if( PRESENT(is_global) )is_global = domain%memory%is_global
    return
  end subroutine mpp_get_memory_domain1d
 
  !#####################################################################
  subroutine mpp_get_compute_domain2d( domain, xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size, &
       x_is_global, y_is_global, tile_count, position )
    type(domain2d),     intent(in) :: domain
    integer, intent(out), optional :: xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size
    logical, intent(out), optional :: x_is_global, y_is_global
    integer, intent(in),  optional :: tile_count, position
    integer                        :: tile, ishift, jshift
 
    tile = 1
    if(present(tile_count)) tile = tile_count
 
    call mpp_get_compute_domain( domain%x(tile), xbegin, xend, xsize, xmax_size, x_is_global )
    call mpp_get_compute_domain( domain%y(tile), ybegin, yend, ysize, ymax_size, y_is_global )
    call mpp_get_domain_shift( domain, ishift, jshift, position )
    if( PRESENT(xend) ) xend  = xend + ishift
    if( PRESENT(yend) ) yend  = yend + jshift
    if( PRESENT(xsize)) xsize = xsize + ishift
    if( PRESENT(ysize)) ysize = ysize + jshift
    if(PRESENT(xmax_size))xmax_size = xmax_size + ishift
    if(PRESENT(ymax_size))ymax_size = ymax_size + jshift
 
    return
  end subroutine mpp_get_compute_domain2d
 
  !#####################################################################
  subroutine mpp_get_data_domain2d( domain, xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size, &
       x_is_global, y_is_global, tile_count, position )
    type(domain2d),     intent(in) :: domain
    integer, intent(out), optional :: xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size
    logical, intent(out), optional :: x_is_global, y_is_global
    integer, intent(in),  optional :: tile_count, position
    integer                        :: tile, ishift, jshift
 
    tile = 1
    if(present(tile_count)) tile = tile_count
 
    call mpp_get_data_domain( domain%x(tile), xbegin, xend, xsize, xmax_size, x_is_global )
    call mpp_get_data_domain( domain%y(tile), ybegin, yend, ysize, ymax_size, y_is_global )
    call mpp_get_domain_shift( domain, ishift, jshift, position )
    if( PRESENT(xend) ) xend  = xend + ishift
    if( PRESENT(yend) ) yend  = yend + jshift
    if( PRESENT(xsize)) xsize = xsize + ishift
    if( PRESENT(ysize)) ysize = ysize + jshift
    if(PRESENT(xmax_size))xmax_size = xmax_size + ishift
    if(PRESENT(ymax_size))ymax_size = ymax_size + jshift
 
    return
  end subroutine mpp_get_data_domain2d
 
  !#####################################################################
  subroutine mpp_get_global_domain2d( domain, xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size, &
                                      tile_count, position )
    type(domain2d),     intent(in) :: domain
    integer, intent(out), optional :: xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size
    integer, intent(in),  optional :: tile_count, position
    integer                        :: tile, ishift, jshift
 
    tile = 1
    if(present(tile_count)) tile = tile_count
 
    call mpp_get_global_domain( domain%x(tile), xbegin, xend, xsize, xmax_size )
    call mpp_get_global_domain( domain%y(tile), ybegin, yend, ysize, ymax_size )
    call mpp_get_domain_shift( domain, ishift, jshift, position )
    if( PRESENT(xend) ) xend  = xend + ishift
    if( PRESENT(yend) ) yend  = yend + jshift
    if( PRESENT(xsize)) xsize = xsize + ishift
    if( PRESENT(ysize)) ysize = ysize + jshift
    if(PRESENT(xmax_size))xmax_size = xmax_size + ishift
    if(PRESENT(ymax_size))ymax_size = ymax_size + jshift
 
    return
  end subroutine mpp_get_global_domain2d
 
  !#####################################################################
  subroutine mpp_get_memory_domain2d( domain, xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size, &
       x_is_global, y_is_global, position)
    type(domain2d),     intent(in) :: domain
    integer, intent(out), optional :: xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size
    logical, intent(out), optional :: x_is_global, y_is_global
    integer, intent(in),  optional :: position
    integer                        :: tile, ishift, jshift
 
    tile = 1
 
    call mpp_get_memory_domain( domain%x(tile), xbegin, xend, xsize, xmax_size, x_is_global )
    call mpp_get_memory_domain( domain%y(tile), ybegin, yend, ysize, ymax_size, y_is_global )
    call mpp_get_domain_shift( domain, ishift, jshift, position )
    if( PRESENT(xend) ) xend  = xend + ishift
    if( PRESENT(yend) ) yend  = yend + jshift
    if( PRESENT(xsize)) xsize = xsize + ishift
    if( PRESENT(ysize)) ysize = ysize + jshift
    if(PRESENT(xmax_size))xmax_size = xmax_size + ishift
    if(PRESENT(ymax_size))ymax_size = ymax_size + jshift
 
    return
  end subroutine mpp_get_memory_domain2d
 
  !> @brief Modifies the indices in the domain_axis_spec type to those of the supergrid
  subroutine mpp_set_super_grid_indices(grid)
   type(domain_axis_spec), intent(inout) :: grid !< domain_axis_spec type
 
   grid%begin = 2*grid%begin-1
   grid%end = 2*grid%end+1
   grid%size = grid%end-grid%begin+1
 
  end subroutine mpp_set_super_grid_indices
 
  !> @brief Modifies the indices of the input domain to create the supergrid domain
  !!
  !> This is an example of how to use mpp_create_super_grid_domain
  !! @code{.F90}
  !! call mpp_copy_domain(domain_in, domain_out)
  !! call super_grid_domain(domain_out)
  !! @endcode
  !! domain_in is the original domain, domain_out is the domain with the supergrid indices.
  subroutine mpp_create_super_grid_domain(domain)
    type(domain2d), intent(inout) :: domain !< Input domain
 
    integer :: xbegin !< Begining x indices
    integer :: ybegin !< Begining y indices
    integer :: xend   !< Ending x indices
    integer :: yend   !< Ending y indices
    integer :: xsize  !< Size of the x domain
    integer :: ysize  !< Size of the y domain
    integer :: i      !< For loops
 
    call mpp_get_global_domain(domain, xbegin=xbegin, xend=xend, ybegin=ybegin, yend=yend, xsize=xsize, ysize=ysize)
    call mpp_set_global_domain (domain, 2*xbegin-1, 2*xend+1, 2*ybegin-1, 2*yend+1, 2*(xend-xbegin)+3, &
                               &  2*(yend-ybegin)+3)
 
    call mpp_get_compute_domain(domain, xbegin=xbegin, xend=xend, ybegin=ybegin, yend=yend, xsize=xsize, ysize=ysize)
    call mpp_set_compute_domain (domain, 2*xbegin-1, 2*xend+1, 2*ybegin-1, 2*yend+1, 2*(xend-xbegin)+3, &
                                &  2*(yend-ybegin)+3)
 
    call mpp_get_data_domain(domain, xbegin=xbegin, xend=xend, ybegin=ybegin, yend=yend, xsize=xsize, ysize=ysize)
    call mpp_set_data_domain (domain, 2*xbegin-1, 2*xend+1, 2*ybegin-1, 2*yend+1, 2*(xend-xbegin)+3, 2*(yend-ybegin)+3)
 
    do i=1, size(domain%list(:))
        call mpp_set_super_grid_indices(domain%list(i-1)%x(1)%global)
        call mpp_set_super_grid_indices(domain%list(i-1)%y(1)%global)
 
        call mpp_set_super_grid_indices(domain%list(i-1)%x(1)%compute)
        call mpp_set_super_grid_indices(domain%list(i-1)%y(1)%compute)
 
        !< There is no data domain in domain%list
        !call mpp_set_super_grid_indices(domain%list(i-1)%x(1)%domain_data)
        !call mpp_set_super_grid_indices(domain%list(i-1)%y(1)%domain_data)
    enddo
 
    do i=1, size(domain%x(1)%list)
        call mpp_set_super_grid_indices(domain%x(1)%list(i-1)%compute)
    enddo
 
    do i=1, size(domain%y(1)%list)
        call mpp_set_super_grid_indices(domain%y(1)%list(i-1)%compute)
    enddo
 
  end subroutine mpp_create_super_grid_domain
 
  !#####################################################################
  subroutine mpp_set_compute_domain1d( domain, begin, end, size, is_global )
    type(domain1d), intent(inout) :: domain
    integer, intent(in), optional :: begin, end, size
    logical, intent(in), optional :: is_global
 
    if(present(begin)) domain%compute%begin = begin
    if(present(end))   domain%compute%end   = end
    if(present(size))  domain%compute%size  = size
    if(present(is_global)) domain%compute%is_global = is_global
 
  end subroutine mpp_set_compute_domain1d
 
  !#####################################################################
  subroutine mpp_set_compute_domain2d( domain, xbegin, xend, ybegin, yend, xsize, ysize, &
                                       x_is_global, y_is_global, tile_count )
    type(domain2d), intent(inout) :: domain
    integer, intent(in), optional :: xbegin, xend, ybegin, yend, xsize, ysize
    logical, intent(in), optional :: x_is_global, y_is_global
    integer, intent(in),  optional :: tile_count
    integer                        :: tile
 
    tile = 1
    if(present(tile_count)) tile = tile_count
 
    call mpp_set_compute_domain(domain%x(tile), xbegin, xend, xsize, x_is_global)
    call mpp_set_compute_domain(domain%y(tile), ybegin, yend, ysize, y_is_global)
 
  end subroutine mpp_set_compute_domain2d
 
  !#####################################################################
  subroutine mpp_set_data_domain1d( domain, begin, end, size, is_global )
    type(domain1d), intent(inout) :: domain
    integer, intent(in), optional :: begin, end, size
    logical, intent(in), optional :: is_global
 
    if(present(begin)) domain%domain_data%begin = begin
    if(present(end))   domain%domain_data%end   = end
    if(present(size))  domain%domain_data%size  = size
    if(present(is_global)) domain%domain_data%is_global = is_global
 
  end subroutine mpp_set_data_domain1d
 
  !#####################################################################
  subroutine mpp_set_data_domain2d( domain, xbegin, xend, ybegin, yend, xsize, ysize, &
                                    x_is_global, y_is_global, tile_count )
    type(domain2d), intent(inout) :: domain
    integer, intent(in), optional :: xbegin, xend, ybegin, yend, xsize, ysize
    logical, intent(in), optional :: x_is_global, y_is_global
    integer, intent(in),  optional :: tile_count
    integer                        :: tile
 
    tile = 1
    if(present(tile_count)) tile = tile_count
 
    call mpp_set_data_domain(domain%x(tile), xbegin, xend, xsize, x_is_global)
    call mpp_set_data_domain(domain%y(tile), ybegin, yend, ysize, y_is_global)
 
  end subroutine mpp_set_data_domain2d
 
  !#####################################################################
  subroutine mpp_set_global_domain1d( domain, begin, end, size)
    type(domain1d), intent(inout) :: domain
    integer, intent(in), optional :: begin, end, size
 
    if(present(begin)) domain%global%begin = begin
    if(present(end))   domain%global%end   = end
    if(present(size))  domain%global%size  = size
 
  end subroutine mpp_set_global_domain1d
 
  !#####################################################################
  subroutine mpp_set_global_domain2d( domain, xbegin, xend, ybegin, yend, xsize, ysize, tile_count )
    type(domain2d), intent(inout) :: domain
    integer, intent(in), optional :: xbegin, xend, ybegin, yend, xsize, ysize
    integer, intent(in),  optional :: tile_count
    integer                        :: tile
 
    tile = 1
    if(present(tile_count)) tile = tile_count
    call mpp_set_global_domain(domain%x(tile), xbegin, xend, xsize)
    call mpp_set_global_domain(domain%y(tile), ybegin, yend, ysize)
 
  end subroutine mpp_set_global_domain2d
 
  !#####################################################################
 
  !> @brief Retrieve 1D components of 2D decomposition.
  !!
  !! It is sometime necessary to have direct recourse to the domain1D types
  !! that compose a domain2D object. This call retrieves them.
  !! @code{.F90}
  !!    call mpp_get_domain_components( domain, x, y )
  !! @endcode
  subroutine mpp_get_domain_components( domain, x, y, tile_count )
    type(domain2d),            intent(in) :: domain
    type(domain1d), intent(inout), optional :: x, y
    integer, intent(in),  optional :: tile_count
    integer                        :: tile
 
    tile = 1
    if(present(tile_count)) tile = tile_count
    if( PRESENT(x) )x = domain%x(tile)
    if( PRESENT(y) )y = domain%y(tile)
    return
  end subroutine mpp_get_domain_components
 
  !#####################################################################
  subroutine mpp_get_compute_domains1d( domain, begin, end, size )
    type(domain1d),                   intent(in) :: domain
    integer, intent(out), optional, dimension(:) :: begin, end, size
 
    if( .NOT.module_is_initialized ) &
         call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: must first call mpp_domains_init.' )
    !we use shape instead of size for error checks because size is used as an argument
    if( PRESENT(begin) )then
       if( any(shape(begin).NE.shape(domain%list)) ) &
            call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: begin array size does not match domain.' )
       begin(:) = domain%list(:)%compute%begin
    end if
    if( PRESENT(end) )then
       if( any(shape(end).NE.shape(domain%list)) ) &
            call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: end array size does not match domain.' )
            end(:) = domain%list(:)%compute%end
    end if
    if( PRESENT(size) )then
       if( any(shape(size).NE.shape(domain%list)) ) &
           call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: size array size does not match domain.' )
           size(:) = domain%list(:)%compute%size
    end if
    return
end subroutine mpp_get_compute_domains1d
 
!#####################################################################
subroutine mpp_get_compute_domains2d( domain, xbegin, xend, xsize, ybegin, yend, ysize, position )
 type(domain2d),                   intent(in) :: domain
 integer, intent(out), optional, dimension(:) :: xbegin, xend, xsize, ybegin, yend, ysize
 integer, intent(in ), optional               :: position
 
 integer :: i, ishift, jshift
 
 call mpp_get_domain_shift( domain, ishift, jshift, position )
 
 
 if( .NOT.module_is_initialized ) &
      call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: must first call mpp_domains_init.' )
 
 if( PRESENT(xbegin) )then
    if( size(xbegin(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: xbegin array size does not match domain.' )
    do i = 1, size(xbegin(:))
       xbegin(i) = domain%list(i-1)%x(1)%compute%begin
    end do
 end if
 if( PRESENT(xend) )then
    if( size(xend(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: xend array size does not match domain.' )
    do i = 1, size(xend(:))
       xend(i) = domain%list(i-1)%x(1)%compute%end + ishift
    end do
 end if
 if( PRESENT(xsize) )then
    if( size(xsize(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: xsize array size does not match domain.' )
    do i = 1, size(xsize(:))
       xsize(i) = domain%list(i-1)%x(1)%compute%size + ishift
    end do
 end if
 if( PRESENT(ybegin) )then
    if( size(ybegin(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: ybegin array size does not match domain.' )
    do i = 1, size(ybegin(:))
       ybegin(i) = domain%list(i-1)%y(1)%compute%begin
    end do
 end if
 if( PRESENT(yend) )then
    if( size(yend(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: yend array size does not match domain.' )
    do i = 1, size(yend(:))
       yend(i) = domain%list(i-1)%y(1)%compute%end + jshift
    end do
 end if
 if( PRESENT(ysize) )then
    if( size(ysize(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_COMPUTE_DOMAINS: ysize array size does not match domain.' )
    do i = 1, size(ysize(:))
       ysize(i) = domain%list(i-1)%y(1)%compute%size + jshift
    end do
 end if
 return
end subroutine mpp_get_compute_domains2d
 
  !#####################################################################
  subroutine mpp_get_global_domains1d( domain, begin, end, size )
    type(domain1d),                   intent(in) :: domain
    integer, intent(out), optional, dimension(:) :: begin, end, size
 
    if( .NOT.module_is_initialized ) &
         call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: must first call mpp_domains_init.' )
    !we use shape instead of size for error checks because size is used as an argument
    if( PRESENT(begin) )then
       if( any(shape(begin).NE.shape(domain%list)) ) &
            call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: begin array size does not match domain.' )
       begin(:) = domain%list(:)%global%begin
    end if
    if( PRESENT(end) )then
       if( any(shape(end).NE.shape(domain%list)) ) &
            call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: end array size does not match domain.' )
            end(:) = domain%list(:)%global%end
    end if
    if( PRESENT(size) )then
       if( any(shape(size).NE.shape(domain%list)) ) &
           call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: size array size does not match domain.' )
           size(:) = domain%list(:)%global%size
    end if
    return
  end subroutine mpp_get_global_domains1d
 
 
!#####################################################################
subroutine mpp_get_global_domains2d( domain, xbegin, xend, xsize, ybegin, yend, ysize, position )
 type(domain2d),                   intent(in) :: domain
 integer, intent(out), optional, dimension(:) :: xbegin, xend, xsize, ybegin, yend, ysize
 integer, intent(in ), optional               :: position
 
 integer :: i, ishift, jshift
 
 call mpp_get_domain_shift( domain, ishift, jshift, position )
 
 
 if( .NOT.module_is_initialized ) &
      call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: must first call mpp_domains_init.' )
 
 if( PRESENT(xbegin) )then
    if( size(xbegin(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: xbegin array size does not match domain.' )
    do i = 1, size(xbegin(:))
       xbegin(i) = domain%list(i-1)%x(1)%global%begin
    end do
 end if
 if( PRESENT(xend) )then
    if( size(xend(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: xend array size does not match domain.' )
    do i = 1, size(xend(:))
       xend(i) = domain%list(i-1)%x(1)%global%end + ishift
    end do
 end if
 if( PRESENT(xsize) )then
    if( size(xsize(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: xsize array size does not match domain.' )
    do i = 1, size(xsize(:))
       xsize(i) = domain%list(i-1)%x(1)%global%size + ishift
    end do
 end if
 if( PRESENT(ybegin) )then
    if( size(ybegin(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: ybegin array size does not match domain.' )
    do i = 1, size(ybegin(:))
       ybegin(i) = domain%list(i-1)%y(1)%global%begin
    end do
 end if
 if( PRESENT(yend) )then
    if( size(yend(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: yend array size does not match domain.' )
    do i = 1, size(yend(:))
       yend(i) = domain%list(i-1)%y(1)%global%end + jshift
    end do
 end if
 if( PRESENT(ysize) )then
    if( size(ysize(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_GLOBAL_DOMAINS: ysize array size does not match domain.' )
    do i = 1, size(ysize(:))
       ysize(i) = domain%list(i-1)%y(1)%global%size + jshift
    end do
 end if
 return
end subroutine mpp_get_global_domains2d
 
 
!#####################################################################
subroutine mpp_get_domain_extents1d(domain, xextent, yextent)
   type(domain2d),            intent(in) :: domain
   integer, dimension(0:), intent(inout) :: xextent, yextent
   integer                               :: n
 
   if(domain%ntiles .NE. 1) call mpp_error(fatal,"mpp_domains_util.inc(mpp_get_domain_extents1D): "// &
       "ntiles is more than 1, please use mpp_get_domain_extents2D")
   if(size(xextent) .NE. size(domain%x(1)%list(:)))  call mpp_error(fatal, &
       & "mpp_domains_util.inc(mpp_get_domain_extents1D): "// &
       & "size(xextent) does not equal to size(domain%x(1)%list(:)))")
   if(size(yextent) .NE. size(domain%y(1)%list(:)))  call mpp_error(fatal, &
       & "mpp_domains_util.inc(mpp_get_domain_extents1D): "// &
       & "size(yextent) does not equal to size(domain%y(1)%list(:)))")
   do n = 0, size(domain%x(1)%list(:))-1
      xextent(n) = domain%x(1)%list(n)%compute%size
   enddo
   do n = 0, size(domain%y(1)%list(:))-1
      yextent(n) = domain%y(1)%list(n)%compute%size
   enddo
 
end subroutine mpp_get_domain_extents1d
 
!#####################################################################
!>  This will return xextent and yextent for each tile
subroutine mpp_get_domain_extents2d(domain, xextent, yextent)
   type(domain2d),             intent(in) :: domain
   integer, dimension(:,:), intent(inout) :: xextent, yextent
   integer                                :: ntile, nlist, n, m, ndivx, ndivy, tile, pos
 
   ntile = domain%ntiles
   nlist = size(domain%list(:))
   if(size(xextent,2) .ne. ntile .or. size(yextent,2) .ne. ntile) call mpp_error(fatal, &
       "mpp_domains_utile.inc: the second dimension size of xextent/yextent is not correct")
   ndivx = size(xextent,1); ndivy = size(yextent,1)
   do n = 0, nlist-1
      if(any(domain%list(n)%x(:)%pos>ndivx-1) ) call mpp_error(fatal, &
         "mpp_domains_utile.inc: first dimension size of xextent is less than the x-layout in some tile")
      if(any(domain%list(n)%y(:)%pos>ndivy-1) ) call mpp_error(fatal, &
         "mpp_domains_utile.inc: first dimension size of yextent is less than the y-layout in some tile")
   end do
 
   xextent = 0; yextent=0
 
   do n = 0, nlist-1
      do m = 1, size(domain%list(n)%tile_id(:))
         tile = domain%list(n)%tile_id(m)
         pos = domain%list(n)%x(m)%pos+1
         if(xextent(pos, tile) == 0) xextent(pos,tile) = domain%list(n)%x(m)%compute%size
         pos = domain%list(n)%y(m)%pos+1
         if(yextent(pos, tile) == 0) yextent(pos,tile) = domain%list(n)%y(m)%compute%size
      end do
   end do
 
 
end subroutine mpp_get_domain_extents2d
 
!#####################################################################
function mpp_get_domain_pe(domain)
   type(domain2d), intent(in) :: domain
   integer                    :: mpp_get_domain_pe
 
   mpp_get_domain_pe = domain%pe
 
 
end function mpp_get_domain_pe
 
 
function mpp_get_domain_tile_root_pe(domain)
   type(domain2d), intent(in) :: domain
   integer                    :: mpp_get_domain_tile_root_pe
 
   mpp_get_domain_tile_root_pe = domain%tile_root_pe
 
end function mpp_get_domain_tile_root_pe
 
 
function mpp_get_domain_tile_commid(domain)
   type(domain2d), intent(in) :: domain                     !> domain you are querying for information
   integer                    :: mpp_get_domain_tile_commid !> declaration of the return tile communicator
 
   mpp_get_domain_tile_commid = domain%tile_comm_id
 
end function mpp_get_domain_tile_commid
 
 
function mpp_get_domain_commid(domain)
   type(domain2d), intent(in) :: domain                !> domain you are querying for information
   integer                    :: mpp_get_domain_commid !> declaration of the return domain communicator
 
   mpp_get_domain_commid = domain%comm_id
 
end function mpp_get_domain_commid
 
 
function mpp_get_io_domain(domain)
   type(domain2d), intent(in) :: domain
   type(domain2d), pointer    :: mpp_get_io_domain
 
   if(ASSOCIATED(domain%io_domain)) then
      mpp_get_io_domain => domain%io_domain
   else
      mpp_get_io_domain => null()
   endif
 
end function mpp_get_io_domain
 
!#####################################################################
! <SUBROUTINE NAME="mpp_get_pelist1D" INTERFACE="mpp_get_pelist">
!  <IN NAME="domain" TYPE="type(domain1D)"></IN>
!  <OUT NAME="pelist" TYPE="integer" DIM="(:)"></OUT>
!  <OUT NAME="pos" TYPE="integer"></OUT>
! </SUBROUTINE>
subroutine mpp_get_pelist1d( domain, pelist, pos )
 type(domain1d),     intent(in) :: domain
 integer,           intent(out) :: pelist(:)
 integer, intent(out), optional :: pos
 integer                        :: ndivs
 
 if( .NOT.module_is_initialized ) &
      call mpp_error( fatal, 'MPP_GET_PELIST: must first call mpp_domains_init.' )
 ndivs = size(domain%list(:))
 
 if( size(pelist(:)).NE.ndivs ) &
      call mpp_error( fatal, 'MPP_GET_PELIST: pelist array size does not match domain.' )
 
 pelist(:) = domain%list(0:ndivs-1)%pe
 if( PRESENT(pos) )pos = domain%pos
 return
end subroutine mpp_get_pelist1d
 
!#####################################################################
! <SUBROUTINE NAME="mpp_get_pelist2D" INTERFACE="mpp_get_pelist">
!  <IN NAME="domain" TYPE="type(domain2D)"></IN>
!  <OUT NAME="pelist" TYPE="integer" DIM="(:)"></OUT>
!  <OUT NAME="pos" TYPE="integer"></OUT>
! </SUBROUTINE>
subroutine mpp_get_pelist2d( domain, pelist, pos )
 type(domain2d),     intent(in) :: domain
 integer,           intent(out) :: pelist(:)
 integer, intent(out), optional :: pos
 
 if( .NOT.module_is_initialized ) &
      call mpp_error( fatal, 'MPP_GET_PELIST: must first call mpp_domains_init.' )
 if( size(pelist(:)).NE.size(domain%list(:)) ) &
      call mpp_error( fatal, 'MPP_GET_PELIST: pelist array size does not match domain.' )
 
 pelist(:) = domain%list(:)%pe
 if( PRESENT(pos) )pos = domain%pos
 return
end subroutine mpp_get_pelist2d
 
!#####################################################################
! <SUBROUTINE NAME="mpp_get_layout1D" INTERFACE="mpp_get_layout">
!  <IN NAME="domain" TYPE="type(domain1D)"></IN>
!  <OUT NAME="layout" TYPE="integer"></OUT>
! </SUBROUTINE>
subroutine mpp_get_layout1d( domain, layout )
 type(domain1d), intent(in) :: domain
 integer,       intent(out) :: layout
 
 if( .NOT.module_is_initialized ) &
      call mpp_error( fatal, 'MPP_GET_LAYOUT: must first call mpp_domains_init.' )
 
 layout = size(domain%list(:))
 return
end subroutine mpp_get_layout1d
 
!#####################################################################
! <SUBROUTINE NAME="mpp_get_layout2D" INTERFACE="mpp_get_layout">
!  <IN NAME="domain" TYPE="type(domain2D)"></IN>
!  <OUT NAME="layout" TYPE="integer" DIM="(2)"></OUT>
! </SUBROUTINE>
subroutine mpp_get_layout2d( domain, layout )
 type(domain2d), intent(in) :: domain
 integer,       intent(out) :: layout(2)
 
 if( .NOT.module_is_initialized ) &
      call mpp_error( fatal, 'MPP_GET_LAYOUT: must first call mpp_domains_init.' )
 
 layout(1) = size(domain%x(1)%list(:))
 layout(2) = size(domain%y(1)%list(:))
 return
end subroutine mpp_get_layout2d
 
!#####################################################################
 
  !> @brief Returns the shift value in x and y-direction according to domain position..
  !!
  !> When domain is symmetry, one extra point maybe needed in
  !!    x- and/or y-direction. This routine will return the shift value based
  !!    on the position
  !! @code{.F90}
  !!   call mpp_get_domain_shift( domain, ishift, jshift, position )
  !! @endcode
subroutine mpp_get_domain_shift(domain, ishift, jshift, position)
  type(domain2d),     intent(in) :: domain!> predefined data contains 2-d domain decomposition.
  integer,           intent(out) :: ishift, jshift!< return value will be 0 or 1.
  integer, optional,  intent(in) :: position!< position of data. Its value can be CENTER, EAST, NORTH or CORNER.
  integer :: pos
 
  ishift = 0 ; jshift = 0
  pos = center
  if(present(position)) pos = position
 
  if(domain%symmetry) then ! shift is non-zero only when the domain is symmetry.
     select case(pos)
     case(corner)
        ishift = 1; jshift = 1
     case(east)
        ishift = 1
     case(north)
        jshift = 1
     end select
  end if
 
end subroutine mpp_get_domain_shift
 
!#####################################################################
 
    !> Return PE to the righ/left of this PE-domain.
    subroutine mpp_get_neighbor_pe_1d(domain, direction, pe)
      type(domain1d), intent(inout) :: domain
      integer, intent(in)           :: direction
      integer, intent(out)          :: pe
 
      integer ipos, ipos2, npx
 
      pe   = null_pe
      npx  = size(domain%list(:)) ! 0..npx-1
      ipos = domain%pos
 
      select case (direction)
 
      case (:-1)
         ! neighbor on the left
         ipos2 = ipos - 1
         if(ipos2 <    0) then
            if(domain%cyclic) then
               ipos2 = npx-1
            else
               ipos2 = -999
            endif
         endif
 
      case (0)
         ! identity
         ipos2 = ipos
 
      case (1:)
         ! neighbor on the right
         ipos2 = ipos + 1
         if(ipos2 > npx-1) then
            if(domain%cyclic) then
               ipos2 = 0
            else
               ipos2 = -999
            endif
         endif
 
      end select
 
      if(ipos2 >= 0) pe = domain%list(ipos2)%pe
 
    end subroutine mpp_get_neighbor_pe_1d
!#####################################################################
 
    !> Return PE North/South/East/West of this PE-domain.
    !! direction must be NORTH, SOUTH, EAST or WEST.
    subroutine mpp_get_neighbor_pe_2d(domain, direction, pe)
      type(domain2d), intent(inout) :: domain
      integer, intent(in)           :: direction
      integer, intent(out)          :: pe
 
      integer ipos, jpos, npx, npy, ix, iy, ipos0, jpos0
 
      pe   = null_pe
      npx  = size(domain%x(1)%list(:)) ! 0..npx-1
      npy  = size(domain%y(1)%list(:)) ! 0..npy-1
      ipos0 = domain%x(1)%pos
      jpos0 = domain%y(1)%pos
 
      select case (direction)
      case (north)
         ix = 0
         iy = 1
      case (north_east)
         ix = 1
         iy = 1
      case (east)
         ix = 1
         iy = 0
      case (south_east)
         ix = 1
         iy =-1
      case (south)
         ix = 0
         iy =-1
      case (south_west)
         ix =-1
         iy =-1
      case (west)
         ix =-1
         iy = 0
      case (north_west)
         ix =-1
         iy = 1
 
      case default
         call mpp_error( fatal, &
 & 'MPP_GET_NEIGHBOR_PE_2D: direction must be either NORTH, ' &
 & // 'SOUTH, EAST, WEST, NORTH_EAST, SOUTH_EAST, SOUTH_WEST or NORTH_WEST')
      end select
 
      ipos = ipos0 + ix
      jpos = jpos0 + iy
 
 
      if( (ipos < 0 .or. ipos > npx-1) .and. domain%x(1)%cyclic ) then
         ! E/W cyclic domain
         ipos = modulo(ipos, npx)
      endif
 
      if(    (ipos < 0     .and. btest(domain%fold,west)) .or. &
           & (ipos > npx-1 .and. btest(domain%fold,east)) ) then
         ! E or W folded domain
           ipos = ipos0
           jpos = npy-jpos-1
      endif
 
      if( (jpos < 0 .or. jpos > npy-1) .and. domain%y(1)%cyclic ) then
         ! N/S cyclic
         jpos = modulo(jpos, npy)
      endif
 
      if(    (jpos < 0     .and. btest(domain%fold,south)) .or. &
           & (jpos > npy-1 .and. btest(domain%fold,north)) ) then
         ! N or S folded
           ipos = npx-ipos-1
           jpos = jpos0
      endif
 
      ! get the PE number
      pe = null_pe
      if(ipos >= 0 .and. ipos <= npx-1 .and. jpos >= 0 .and. jpos <= npy-1) then
         pe = domain%pearray(ipos, jpos)
      endif
 
 
    end subroutine mpp_get_neighbor_pe_2d
 
 
!#######################################################################
 
  subroutine nullify_domain2d_list(domain)
     type(domain2d), intent(inout) :: domain
 
     domain%list =>null()
 
  end subroutine nullify_domain2d_list
 
!#######################################################################
  function mpp_domain_is_symmetry(domain)
    type(domain2d), intent(in) :: domain
    logical                    :: mpp_domain_is_symmetry
 
    mpp_domain_is_symmetry = domain%symmetry
    return
 
  end function mpp_domain_is_symmetry
 
!#######################################################################
  function mpp_domain_is_initialized(domain)
    type(domain2d), intent(in) :: domain
    logical                    :: mpp_domain_is_initialized
 
    mpp_domain_is_initialized = domain%initialized
 
    return
 
  end function mpp_domain_is_initialized
 
!#######################################################################
  !--- private routine used only for mpp_update_domains. This routine will
  !--- compare whalo, ehalo, shalo, nhalo with the halo size when defining "domain"
  !--- to decide if update is needed. Also it check the sign of whalo, ehalo, shalo and nhalo.
  function domain_update_is_needed(domain, whalo, ehalo, shalo, nhalo)
    type(domain2d), intent(in) :: domain
    integer,        intent(in) :: whalo, ehalo, shalo, nhalo
    logical                    :: domain_update_is_needed
 
    domain_update_is_needed = .true.
 
    if(whalo == 0 .AND. ehalo==0 .AND. shalo == 0 .AND. nhalo==0 ) then
       domain_update_is_needed = .false.
       if( debug )call mpp_error(note, &
          'mpp_domains_util.inc: halo size to be updated are all zero, no update will be done')
       return
    end if
    if(  (whalo == -domain%whalo .AND. domain%whalo .NE. 0) .or.  &
         (ehalo == -domain%ehalo .AND. domain%ehalo .NE. 0) .or.  &
         (shalo == -domain%shalo .AND. domain%shalo .NE. 0) .or.  &
         (nhalo == -domain%nhalo .AND. domain%nhalo .NE. 0) ) then
       domain_update_is_needed = .false.
       call mpp_error(note, 'mpp_domains_util.inc: at least one of w/e/s/n halo size to be updated '// &
            'is the inverse of the original halo when defining domain, no update will be done')
       return
    end if
 
  end function domain_update_is_needed
!#######################################################################
  !> this routine found the domain has the same halo size with the input
  !! whalo, ehalo,
  function search_update_overlap(domain, whalo, ehalo, shalo, nhalo, position)
    type(domain2d),          intent(inout) :: domain
    integer,                    intent(in) :: whalo, ehalo, shalo, nhalo
    integer,                    intent(in) :: position
    type(overlapspec),          pointer    :: search_update_overlap
    type(overlapspec),          pointer    :: update_ref
    type(overlapspec), pointer :: check => null()
    integer :: ishift, jshift, shift
 
    shift = 0; if(domain%symmetry) shift = 1
    select case(position)
    case (center)
       update_ref => domain%update_T
       ishift = 0; jshift = 0
    case (corner)
       update_ref => domain%update_C
       ishift = shift; jshift = shift
    case (north)
       update_ref => domain%update_N
       ishift = 0; jshift = shift
    case (east)
       update_ref => domain%update_E
       ishift = shift; jshift = 0
    case default
       call mpp_error(fatal,"mpp_domains_util.inc(search_update_overlap): position should be CENTER|CORNER|EAST|NORTH")
    end select
 
    search_update_overlap => update_ref
 
    do
       if(whalo == search_update_overlap%whalo .AND. ehalo == search_update_overlap%ehalo .AND. &
            shalo == search_update_overlap%shalo .AND. nhalo == search_update_overlap%nhalo ) then
            exit ! found domain
       endif
       !--- if not found, switch to next
       if(.NOT. ASSOCIATED(search_update_overlap%next)) then
          allocate(search_update_overlap%next)
          search_update_overlap => search_update_overlap%next
          if(domain%fold .NE. 0) then
             call compute_overlaps(domain, position, search_update_overlap, check, &
                                   ishift, jshift, 0, 0, whalo, ehalo, shalo, nhalo)
          else
             call set_overlaps(domain, update_ref, search_update_overlap, whalo, ehalo, shalo, nhalo )
          endif
          exit
       else
          search_update_overlap => search_update_overlap%next
       end if
 
    end do
 
    update_ref => null()
 
  end function search_update_overlap
 
!#######################################################################
  !> this routine finds the check at certain position
  function search_check_overlap(domain, position)
    type(domain2d),             intent(in) :: domain
    integer,                    intent(in) :: position
    type(overlapspec),          pointer    :: search_check_overlap
 
    select case(position)
    case (center)
       search_check_overlap => null()
    case (corner)
       search_check_overlap => domain%check_C
    case (north)
       search_check_overlap => domain%check_N
    case (east)
       search_check_overlap => domain%check_E
    case default
       call mpp_error(fatal,"mpp_domains_util.inc(search_check_overlap): position should be CENTER|CORNER|EAST|NORTH")
    end select
 
  end function search_check_overlap
 
!#######################################################################
  !> This routine finds the bound at certain position
  function search_bound_overlap(domain, position)
    type(domain2d),             intent(in) :: domain
    integer,                    intent(in) :: position
    type(overlapspec),          pointer    :: search_bound_overlap
 
    select case(position)
    case (center)
       search_bound_overlap => null()
    case (corner)
       search_bound_overlap => domain%bound_C
    case (north)
       search_bound_overlap => domain%bound_N
    case (east)
       search_bound_overlap => domain%bound_E
    case default
       call mpp_error(fatal,"mpp_domains_util.inc(search_bound_overlap): position should be CENTER|CORNER|EAST|NORTH")
    end select
 
  end function search_bound_overlap
 
  !########################################################################
  !> Returns the tile_id on current pe
  function mpp_get_tile_id(domain)
     type(domain2d),                  intent(in) :: domain
     integer, dimension(size(domain%tile_id(:))) :: mpp_get_tile_id
 
     mpp_get_tile_id = domain%tile_id
     return
 
  end function mpp_get_tile_id
 
  !#######################################################################
  !> Return the tile_id on current pelist. one-tile-per-pe is assumed.
  subroutine mpp_get_tile_list(domain, tiles)
     type(domain2d), intent(in) :: domain
     integer,     intent(inout) :: tiles(:)
     integer                    :: i
 
    if( size(tiles(:)).NE.size(domain%list(:)) ) &
         call mpp_error( fatal, 'mpp_get_tile_list: tiles array size does not match domain.' )
    do i = 1, size(tiles(:))
       if(size(domain%list(i-1)%tile_id(:)) > 1) call mpp_error( fatal,  &
               'mpp_get_tile_list: only support one-tile-per-pe now, contact developer');
       tiles(i) = domain%list(i-1)%tile_id(1)
    end do
 
  end subroutine mpp_get_tile_list
 
  !########################################################################
  !> Returns number of tiles in mosaic
  function mpp_get_ntile_count(domain)
     type(domain2d),                  intent(in) :: domain
     integer                                     :: mpp_get_ntile_count
 
     mpp_get_ntile_count = domain%ntiles
     return
 
  end function mpp_get_ntile_count
 
  !########################################################################
  !> Returns number of tile on current pe
  function mpp_get_current_ntile(domain)
     type(domain2d),                  intent(in) :: domain
     integer                                     :: mpp_get_current_ntile
 
     mpp_get_current_ntile = size(domain%tile_id(:))
     return
 
  end function mpp_get_current_ntile
 
  !#######################################################################
  !> Returns if current pe is the root pe of the tile, if number of tiles on current pe
  !! is greater than 1, will return true, if isc==isg and jsc==jsg also will return true,
  !! otherwise false will be returned.
  function mpp_domain_is_tile_root_pe(domain)
     type(domain2d), intent(in) :: domain
     logical                    :: mpp_domain_is_tile_root_pe
 
     mpp_domain_is_tile_root_pe = domain%pe == domain%tile_root_pe;
 
  end function mpp_domain_is_tile_root_pe
 
  !#########################################################################
  !> Returns number of processors used on current tile.
  function mpp_get_tile_npes(domain)
    type(domain2d), intent(in)  :: domain
    integer                     :: mpp_get_tile_npes
    integer                     :: i, tile
 
    !--- When there is more than one tile on this pe, we assume each tile will be
    !--- limited to this pe.
    if(size(domain%tile_id(:)) > 1) then
       mpp_get_tile_npes = 1
    else
      mpp_get_tile_npes = 0
      tile = domain%tile_id(1)
      do i = 0, size(domain%list(:))-1
         if(tile == domain%list(i)%tile_id(1) )  mpp_get_tile_npes =  mpp_get_tile_npes + 1
      end do
    endif
 
  end function mpp_get_tile_npes
 
  !########################################################################
  !> Get the processors list used on current tile.
  subroutine mpp_get_tile_pelist(domain, pelist)
    type(domain2d), intent(in)  :: domain
    integer,      intent(inout) :: pelist(:)
    integer                     :: npes_on_tile
    integer                     :: i, tile, pos
 
    npes_on_tile = mpp_get_tile_npes(domain)
    if(size(pelist(:)) .NE. npes_on_tile) call mpp_error(fatal, &
           "mpp_domains_util.inc(mpp_get_tile_pelist): size(pelist) does not equal npes on current tile")
    tile = domain%tile_id(1)
    pos = 0
    do i = 0, size(domain%list(:))-1
       if(tile == domain%list(i)%tile_id(1)) then
          pos = pos+1
          pelist(pos) = domain%list(i)%pe
       endif
    enddo
 
    return
 
  end subroutine mpp_get_tile_pelist
 
!#####################################################################
subroutine mpp_get_tile_compute_domains( domain, xbegin, xend, ybegin, yend, position )
 type(domain2d),         intent(in) :: domain
 integer, intent(out), dimension(:) :: xbegin, xend, ybegin, yend
 integer, intent(in ), optional     :: position
 
 integer :: i, ishift, jshift
 integer :: npes_on_tile, pos, tile
 
 call mpp_get_domain_shift( domain, ishift, jshift, position )
 
 
 if( .NOT.module_is_initialized ) &
      call mpp_error( fatal, 'mpp_get_compute_domains2D: must first call mpp_domains_init.' )
 
 npes_on_tile = mpp_get_tile_npes(domain)
 if(size(xbegin(:)) .NE. npes_on_tile) call mpp_error(fatal, &
       "mpp_domains_util.inc(mpp_get_compute_domains2D): size(xbegin) does not equal npes on current tile")
 if(size(xend(:)) .NE. npes_on_tile) call mpp_error(fatal, &
       "mpp_domains_util.inc(mpp_get_compute_domains2D): size(xend) does not equal npes on current tile")
 if(size(ybegin(:)) .NE. npes_on_tile) call mpp_error(fatal, &
       "mpp_domains_util.inc(mpp_get_compute_domains2D): size(ybegin) does not equal npes on current tile")
 if(size(yend(:)) .NE. npes_on_tile) call mpp_error(fatal, &
       "mpp_domains_util.inc(mpp_get_compute_domains2D): size(yend) does not equal npes on current tile")
 
 tile = domain%tile_id(1)
 pos = 0
 do i = 0, size(domain%list(:))-1
    if(tile == domain%list(i)%tile_id(1)) then
       pos = pos+1
       xbegin(pos) = domain%list(i)%x(1)%compute%begin
       xend(pos) = domain%list(i)%x(1)%compute%end + ishift
       ybegin(pos) = domain%list(i)%y(1)%compute%begin
       yend(pos) = domain%list(i)%y(1)%compute%end + jshift
    endif
 enddo
 
 return
 
end subroutine mpp_get_tile_compute_domains
 
 
 
  !#############################################################################
  function mpp_get_num_overlap(domain, action, p, position)
    type(domain2d),    intent(in) :: domain
    integer,           intent(in) :: action
    integer,           intent(in) :: p
    integer, optional, intent(in) :: position
    integer                       :: mpp_get_num_overlap
    type(overlapspec), pointer    :: update => null()
    integer                       :: pos
 
    pos = center
    if(present(position)) pos = position
    select case(pos)
    case (center)
       update => domain%update_T
    case (corner)
       update => domain%update_C
    case (east)
       update => domain%update_E
    case (north)
       update => domain%update_N
    case default
       call mpp_error( fatal, "mpp_domains_mod(mpp_get_num_overlap): invalid option of position")
    end select
 
    if(action == event_send) then
       if(p< 1 .OR. p > update%nsend) call mpp_error( fatal, &
                "mpp_domains_mod(mpp_get_num_overlap): p should be between 1 and update%nsend")
       mpp_get_num_overlap = update%send(p)%count
    else if(action == event_recv) then
       if(p< 1 .OR. p > update%nrecv) call mpp_error( fatal, &
                "mpp_domains_mod(mpp_get_num_overlap): p should be between 1 and update%nrecv")
       mpp_get_num_overlap = update%recv(p)%count
    else
       call mpp_error( fatal, "mpp_domains_mod(mpp_get_num_overlap): invalid option of action")
    end if
 
  end function mpp_get_num_overlap
 
  !#############################################################################
  subroutine mpp_get_update_size(domain, nsend, nrecv, position)
    type(domain2d),      intent(in) :: domain
    integer,            intent(out) :: nsend, nrecv
    integer, optional,   intent(in) :: position
    integer                         :: pos
 
    pos = center
    if(present(position)) pos = position
    select case(pos)
    case (center)
       nsend = domain%update_T%nsend
       nrecv = domain%update_T%nrecv
    case (corner)
       nsend = domain%update_C%nsend
       nrecv = domain%update_C%nrecv
    case (east)
       nsend = domain%update_E%nsend
       nrecv = domain%update_E%nrecv
    case (north)
       nsend = domain%update_N%nsend
       nrecv = domain%update_N%nrecv
    case default
       call mpp_error( fatal, "mpp_domains_mod(mpp_get_update_size): invalid option of position")
    end select
 
  end subroutine mpp_get_update_size
 
  !#############################################################################
  subroutine mpp_get_update_pelist(domain, action, pelist, position)
    type(domain2d),       intent(in) :: domain
    integer,              intent(in) :: action
    integer,           intent(inout) :: pelist(:)
    integer, optional,    intent(in) :: position
    type(overlapspec), pointer    :: update => null()
    integer                       :: pos, p
 
    pos = center
    if(present(position)) pos = position
    select case(pos)
    case (center)
       update => domain%update_T
    case (corner)
       update => domain%update_C
    case (east)
       update => domain%update_E
    case (north)
       update => domain%update_N
    case default
       call mpp_error( fatal, "mpp_domains_mod(mpp_get_update_pelist): invalid option of position")
    end select
 
    if(action == event_send) then
       if(size(pelist) .NE. update%nsend) call mpp_error( fatal,  &
             "mpp_domains_mod(mpp_get_update_pelist): size of pelist does not match update%nsend")
       do p = 1, update%nsend
          pelist(p) = update%send(p)%pe
       enddo
    else if(action == event_recv) then
       if(size(pelist) .NE. update%nrecv) call mpp_error( fatal,  &
             "mpp_domains_mod(mpp_get_update_pelist): size of pelist does not match update%nrecv")
       do p = 1, update%nrecv
          pelist(p) = update%recv(p)%pe
       enddo
    else
       call mpp_error( fatal, "mpp_domains_mod(mpp_get_update_pelist): invalid option of action")
    end if
 
  end subroutine mpp_get_update_pelist
 
  !#############################################################################
  subroutine mpp_get_overlap(domain, action, p, is, ie, js, je, dir, rot, position)
    type(domain2d),         intent(in) :: domain
    integer,                intent(in) :: action
    integer,                intent(in) :: p
    integer, dimension(:), intent(out) :: is, ie, js, je
    integer, dimension(:), intent(out) :: dir, rot
    integer, optional,      intent(in) :: position
    type(overlapspec),      pointer    :: update => null()
    type(overlap_type),     pointer    :: overlap => null()
    integer                            :: count, pos
 
    pos = center
    if(present(position)) pos = position
    select case(pos)
    case (center)
       update => domain%update_T
    case (corner)
       update => domain%update_C
    case (east)
       update => domain%update_E
    case (north)
       update => domain%update_N
    case default
       call mpp_error( fatal, "mpp_domains_mod(mpp_get_overlap): invalid option of position")
    end select
 
    if(action == event_send) then
       overlap => update%send(p)
    else if(action == event_recv) then
       overlap => update%recv(p)
    else
       call mpp_error( fatal, "mpp_domains_mod(mpp_get_overlap): invalid option of action")
    end if
 
    count = overlap%count
    if(size(is(:)) .NE. count .OR. size(ie(:)) .NE. count .OR. size(js(:)) .NE. count .OR.  &
       size(je(:)) .NE. count .OR. size(dir(:)) .NE. count .OR. size(rot(:)) .NE. count ) &
       call mpp_error( fatal, &
                      &  "mpp_domains_mod(mpp_get_overlap): size mismatch between number of overlap and array size")
 
    is  = overlap%is      (1:count)
    ie  = overlap%ie      (1:count)
    js  = overlap%js      (1:count)
    je  = overlap%je      (1:count)
    dir = overlap%dir     (1:count)
    rot = overlap%rotation(1:count)
 
    update  => null()
    overlap => null()
 
  end subroutine mpp_get_overlap
 
  !##################################################################
  function mpp_get_domain_name(domain)
     type(domain2d),    intent(in) :: domain
     character(len=NAME_LENGTH)    :: mpp_get_domain_name
 
     mpp_get_domain_name = domain%name
 
  end function mpp_get_domain_name
 
  !#################################################################
  function mpp_get_domain_root_pe(domain)
     type(domain2d),    intent(in) :: domain
     integer                       :: mpp_get_domain_root_pe
 
     mpp_get_domain_root_pe = domain%list(0)%pe
 
  end function mpp_get_domain_root_pe
  !#################################################################
  function mpp_get_domain_npes(domain)
     type(domain2d),    intent(in) :: domain
     integer                       :: mpp_get_domain_npes
 
     mpp_get_domain_npes = size(domain%list(:))
 
     return
 
  end function mpp_get_domain_npes
 
  !################################################################
  subroutine mpp_get_domain_pelist(domain, pelist)
     type(domain2d),    intent(in) :: domain
     integer,          intent(out) :: pelist(:)
     integer :: p
 
     if(size(pelist(:)) .NE. size(domain%list(:)) ) then
        call mpp_error(fatal, .NE."mpp_get_domain_pelist: size(pelist(:))  size(domain%list(:)) ")
     endif
 
     do p = 0, size(domain%list(:))-1
        pelist(p+1) = domain%list(p)%pe
     enddo
 
     return
 
  end subroutine mpp_get_domain_pelist
 
  !#################################################################
  function mpp_get_io_domain_layout(domain)
     type(domain2d),    intent(in) :: domain
     integer, dimension(2)         :: mpp_get_io_domain_layout
 
     mpp_get_io_domain_layout = domain%io_layout
 
  end function mpp_get_io_domain_layout
 
  !################################################################
  function get_rank_send(domain, overlap_x, overlap_y, rank_x, rank_y, ind_x, ind_y)
    type(domain2d),    intent(in) :: domain
    type(overlapspec), intent(in) :: overlap_x, overlap_y
    integer,          intent(out) :: rank_x, rank_y, ind_x, ind_y
    integer                       :: get_rank_send
    integer                       :: nlist, nsend_x, nsend_y
 
    nlist = size(domain%list(:))
    nsend_x = overlap_x%nsend
    nsend_y = overlap_y%nsend
    rank_x = nlist+1
    rank_y = nlist+1
    if(nsend_x>0) rank_x = overlap_x%send(1)%pe - domain%pe
    if(nsend_y>0) rank_y = overlap_y%send(1)%pe - domain%pe
    if(rank_x .LT. 0) rank_x = rank_x + nlist
    if(rank_y .LT. 0) rank_y = rank_y + nlist
    get_rank_send = min(rank_x, rank_y)
    ind_x = nsend_x + 1
    ind_y = nsend_y + 1
    if(get_rank_send < nlist+1) then
       if(nsend_x>0) ind_x = 1
       if(nsend_y>0) ind_y = 1
    endif
 
  end function get_rank_send
 
  !############################################################################
  function get_rank_recv(domain, overlap_x, overlap_y, rank_x, rank_y, ind_x, ind_y)
    type(domain2d),    intent(in) :: domain
    type(overlapspec), intent(in) :: overlap_x, overlap_y
    integer,          intent(out) :: rank_x, rank_y, ind_x, ind_y
    integer                       :: get_rank_recv
    integer                       :: nlist, nrecv_x, nrecv_y
 
    nlist = size(domain%list(:))
    nrecv_x = overlap_x%nrecv
    nrecv_y = overlap_y%nrecv
    rank_x = -1
    rank_y = -1
    if(nrecv_x>0) then
       rank_x = overlap_x%recv(1)%pe - domain%pe
       if(rank_x .LE. 0) rank_x = rank_x + nlist
    endif
    if(nrecv_y>0) then
       rank_y = overlap_y%recv(1)%pe - domain%pe
       if(rank_y .LE. 0) rank_y = rank_y + nlist
    endif
    get_rank_recv = max(rank_x, rank_y)
    ind_x = nrecv_x + 1
    ind_y = nrecv_y + 1
    if(get_rank_recv < nlist+1) then
       if(nrecv_x>0) ind_x = 1
       if(nrecv_y>0) ind_y = 1
    endif
 
  end function get_rank_recv
 
  function get_vector_recv(domain, update_x, update_y, ind_x, ind_y, start_pos, pelist)
    type(domain2d),    intent(in) :: domain
    type(overlapspec), intent(in) :: update_x, update_y
    integer,          intent(out) :: ind_x(:), ind_y(:)
    integer,          intent(out) :: start_pos(:)
    integer,          intent(out) :: pelist(:)
    integer                       :: nlist, nrecv_x, nrecv_y, ntot, n
    integer                       :: ix, iy, rank_x, rank_y, cur_pos
    integer                       :: get_vector_recv
 
    nlist = size(domain%list(:))
    nrecv_x = update_x%nrecv
    nrecv_y = update_y%nrecv
 
    ntot = nrecv_x + nrecv_y
 
    n  = 1
    ix = 1
    iy = 1
    ind_x = -1
    ind_y = -1
    get_vector_recv = 0
    cur_pos = 0
    do while (n<=ntot)
       if(ix <= nrecv_x ) then
          rank_x = update_x%recv(ix)%pe-domain%pe
          if(rank_x .LE. 0) rank_x = rank_x + nlist
       else
          rank_x = -1
       endif
       if(iy <= nrecv_y ) then
          rank_y = update_y%recv(iy)%pe-domain%pe
          if(rank_y .LE. 0) rank_y = rank_y + nlist
       else
          rank_y = -1
       endif
       get_vector_recv = get_vector_recv + 1
       start_pos(get_vector_recv) = cur_pos
       if( rank_x == rank_y ) then
          n = n+2
          ind_x(get_vector_recv) = ix
          ind_y(get_vector_recv) = iy
          cur_pos = cur_pos + update_x%recv(ix)%totsize + update_y%recv(iy)%totsize
          pelist(get_vector_recv) = update_x%recv(ix)%pe
          ix = ix + 1
          iy = iy + 1
       else if ( rank_x > rank_y ) then
          n = n+1
          ind_x(get_vector_recv) = ix
          ind_y(get_vector_recv) = -1
          cur_pos = cur_pos + update_x%recv(ix)%totsize
          pelist(get_vector_recv) = update_x%recv(ix)%pe
          ix = ix + 1
       else if ( rank_y > rank_x ) then
          n = n+1
          ind_x(get_vector_recv) = -1
          ind_y(get_vector_recv) = iy
          cur_pos = cur_pos + update_y%recv(iy)%totsize
          pelist(get_vector_recv) = update_y%recv(iy)%pe
          iy = iy+1
       endif
    end do
 
 
  end function get_vector_recv
 
  function get_vector_send(domain, update_x, update_y, ind_x, ind_y, start_pos, pelist)
    type(domain2d),    intent(in) :: domain
    type(overlapspec), intent(in) :: update_x, update_y
    integer,          intent(out) :: ind_x(:), ind_y(:)
    integer,          intent(out) :: start_pos(:)
    integer,          intent(out) :: pelist(:)
    integer                       :: nlist, nsend_x, nsend_y, ntot, n
    integer                       :: ix, iy, rank_x, rank_y, cur_pos
    integer                       :: get_vector_send
 
    nlist = size(domain%list(:))
    nsend_x = update_x%nsend
    nsend_y = update_y%nsend
 
    ntot = nsend_x + nsend_y
    n  = 1
    ix = 1
    iy = 1
    ind_x = -1
    ind_y = -1
    get_vector_send = 0
    cur_pos = 0
    do while (n<=ntot)
       if(ix <= nsend_x ) then
          rank_x = update_x%send(ix)%pe-domain%pe
          if(rank_x .LT. 0) rank_x = rank_x + nlist
       else
          rank_x = nlist+1
       endif
       if(iy <= nsend_y ) then
          rank_y = update_y%send(iy)%pe-domain%pe
          if(rank_y .LT. 0) rank_y = rank_y + nlist
       else
          rank_y = nlist+1
       endif
       get_vector_send = get_vector_send + 1
       start_pos(get_vector_send) = cur_pos
 
       if( rank_x == rank_y ) then
          n = n+2
          ind_x(get_vector_send) = ix
          ind_y(get_vector_send) = iy
          cur_pos = cur_pos + update_x%send(ix)%totsize + update_y%send(iy)%totsize
          pelist(get_vector_send) = update_x%send(ix)%pe
          ix = ix + 1
          iy = iy + 1
       else if ( rank_x < rank_y ) then
          n = n+1
          ind_x(get_vector_send) = ix
          ind_y(get_vector_send) = -1
          cur_pos = cur_pos + update_x%send(ix)%totsize
          pelist(get_vector_send) = update_x%send(ix)%pe
          ix = ix + 1
       else if ( rank_y < rank_x ) then
          n = n+1
          ind_x(get_vector_send) = -1
          ind_y(get_vector_send) = iy
          cur_pos = cur_pos + update_y%send(iy)%totsize
          pelist(get_vector_send) = update_y%send(iy)%pe
          iy = iy+1
       endif
    end do
 
 
  end function get_vector_send
 
 
  !############################################################################
  function get_rank_unpack(domain, overlap_x, overlap_y, rank_x, rank_y, ind_x, ind_y)
    type(domain2d),    intent(in) :: domain
    type(overlapspec), intent(in) :: overlap_x, overlap_y
    integer,          intent(out) :: rank_x, rank_y, ind_x, ind_y
    integer                       :: get_rank_unpack
    integer                       :: nlist, nrecv_x, nrecv_y
 
    nlist = size(domain%list(:))
    nrecv_x = overlap_x%nrecv
    nrecv_y = overlap_y%nrecv
 
    rank_x = nlist+1
    rank_y = nlist+1
    if(nrecv_x>0) rank_x = overlap_x%recv(nrecv_x)%pe - domain%pe
    if(nrecv_y>0) rank_y = overlap_y%recv(nrecv_y)%pe - domain%pe
    if(rank_x .LE.0) rank_x = rank_x + nlist
    if(rank_y .LE.0) rank_y = rank_y + nlist
 
    get_rank_unpack = min(rank_x, rank_y)
    ind_x = 0
    ind_y = 0
    if(get_rank_unpack < nlist+1) then
       if(nrecv_x >0) ind_x = nrecv_x
       if(nrecv_y >0) ind_y = nrecv_y
    endif
 
  end function get_rank_unpack
 
  function get_mesgsize(overlap, do_dir)
   type(overlap_type), intent(in) :: overlap
   logical,            intent(in) :: do_dir(:)
   integer                       :: get_mesgsize
   integer                       :: n, dir
 
   get_mesgsize = 0
   do n = 1, overlap%count
      dir = overlap%dir(n)
      if(do_dir(dir)) then
         get_mesgsize = get_mesgsize + overlap%msgsize(n)
      end if
   end do
 
  end function get_mesgsize
 
  !#############################################################################
  subroutine mpp_set_domain_symmetry(domain, symmetry)
    type(domain2d), intent(inout) :: domain
    logical,        intent(in   ) :: symmetry
 
    domain%symmetry = symmetry
 
  end subroutine mpp_set_domain_symmetry
 
  !> @brief Copies input 1d domain to the output 1d domain
  recursive subroutine mpp_copy_domain1d(domain_in, domain_out)
     type(domain1d),    intent(in) :: domain_in  !< Input domain
     type(domain1d), intent(inout) :: domain_out !< Output domain
 
     integer                       :: i        !< For loop
     integer                       :: starting !< Starting bounds
     integer                       :: ending   !< Ending bounds
 
     domain_out%compute = domain_in%compute
     domain_out%domain_data    = domain_in%domain_data
     domain_out%global  = domain_in%global
     domain_out%memory  = domain_in%memory
     domain_out%cyclic  = domain_in%cyclic
     domain_out%pe      = domain_in%pe
     domain_out%pos     = domain_in%pos
 
     if (associated(domain_in%list)) then
         starting = lbound(domain_in%list, 1)
         ending = ubound(domain_in%list, 1)
         if (associated(domain_out%list)) deallocate(domain_out%list) !< Check if allocated
         allocate(domain_out%list(starting:ending))
 
         do i = starting, ending
            call mpp_copy_domain1d(domain_in%list(i), domain_out%list(i))
         enddo
 
     endif
 
  end subroutine mpp_copy_domain1d
 
  !#################################################################
  !> @brief Copies input 2d domain to the output 2d domain
  subroutine mpp_copy_domain2d(domain_in, domain_out)
     type(domain2d),    intent(in) :: domain_in  !< Input domain
     type(domain2d), intent(inout) :: domain_out !< Output domain
 
     integer :: n,i          !< For loops
     integer :: ntiles       !< Number of tiles
     integer :: starting(2)  !< Starting bounds
     integer :: ending(2)    !< Ending bounds
 
     if (associated(domain_out%x)) then
         call mpp_error(fatal, "mpp_copy_domain: domain_out is already set")
     endif
 
     domain_out%id             = domain_in%id
     domain_out%pe             = domain_in%pe
     domain_out%fold           = domain_in%fold
     domain_out%pos            = domain_in%pos
     domain_out%symmetry       = domain_in%symmetry
     domain_out%whalo          = domain_in%whalo
     domain_out%ehalo          = domain_in%ehalo
     domain_out%shalo          = domain_in%shalo
     domain_out%nhalo          = domain_in%nhalo
     domain_out%ntiles         = domain_in%ntiles
     domain_out%max_ntile_pe   = domain_in%max_ntile_pe
     domain_out%ncontacts      = domain_in%ncontacts
     domain_out%rotated_ninety = domain_in%rotated_ninety
     domain_out%initialized    = domain_in%initialized
     domain_out%tile_root_pe   = domain_in%tile_root_pe
     domain_out%io_layout      = domain_in%io_layout
     domain_out%name           = domain_in%name
 
     ntiles = size(domain_in%x(:))
     allocate(domain_out%x(ntiles), domain_out%y(ntiles), domain_out%tile_id(ntiles) )
     do n = 1, ntiles
        call mpp_copy_domain1d(domain_in%x(n), domain_out%x(n))
        call mpp_copy_domain1d(domain_in%y(n), domain_out%y(n))
     enddo
 
     if (associated(domain_in%pearray)) then
         starting = lbound(domain_in%pearray)
         ending = ubound(domain_in%pearray)
 
         allocate(domain_out%pearray(starting(1):ending(1), starting(2):ending(2)))
         domain_out%pearray=domain_in%pearray
     endif
 
     if (associated(domain_in%tile_id)) then
         starting(1) = lbound(domain_in%tile_id,1)
         ending(1) = ubound(domain_in%tile_id,1)
 
         allocate(domain_out%tile_id(starting(1):ending(1)))
         domain_out%tile_id = domain_in%tile_id
     endif
 
     if (associated(domain_in%tile_id_all)) then
         starting(1) = lbound(domain_in%tile_id_all,1)
         ending(1) = ubound(domain_in%tile_id_all,1)
 
         allocate(domain_out%tile_id_all(starting(1):ending(1)))
         domain_out%tile_id_all = domain_in%tile_id_all
     endif
 
     if (associated(domain_in%list)) then
         starting(1) = lbound(domain_in%list,1)
         ending(1) = ubound(domain_in%list,1)
 
         allocate(domain_out%list(starting(1):ending(1)))
         do i = starting(1), ending(1)
            call mpp_copy_domain2d_spec(domain_in%list(i),domain_out%list(i))
         enddo
     endif
 
     return
 
  end subroutine mpp_copy_domain2d
 
  !> @brief Copies input 2d domain spec to the output 2d domain spec
  subroutine mpp_copy_domain2d_spec(domain2D_spec_in, domain2d_spec_out)
  type(domain2d_spec), intent(in) :: domain2D_spec_in   !< Input
  type(domain2d_spec), intent(out) :: domain2D_spec_out !< Output
 
  integer :: starting  !< Starting bounds
  integer :: ending    !< Ending bounds
  integer :: i         !< For loop
 
  domain2d_spec_out%pe = domain2d_spec_in%pe
  domain2d_spec_out%pos = domain2d_spec_in%pos
  domain2d_spec_out%tile_root_pe = domain2d_spec_in%tile_root_pe
 
  if (associated(domain2d_spec_in%tile_id)) then
      starting = lbound(domain2d_spec_in%tile_id,1)
      ending = ubound(domain2d_spec_in%tile_id,1)
 
      if (associated(domain2d_spec_out%tile_id)) deallocate(domain2d_spec_out%tile_id) !< Check if allocated
      allocate(domain2d_spec_out%tile_id(starting:ending))
      domain2d_spec_out%tile_id = domain2d_spec_in%tile_id
  endif
 
  if (associated(domain2d_spec_in%x)) then
      starting = lbound(domain2d_spec_in%x,1)
      ending = ubound(domain2d_spec_in%x,1)
 
      if (associated(domain2d_spec_out%x)) deallocate(domain2d_spec_out%x) !< Check if allocated
      allocate(domain2d_spec_out%x(starting:ending))
      do i = starting, ending
         call mpp_copy_domain1d_spec(domain2d_spec_in%x(i), domain2d_spec_out%x(i))
      enddo
  endif
 
  if (associated(domain2d_spec_in%y)) then
      starting = lbound(domain2d_spec_in%y,1)
      ending = ubound(domain2d_spec_in%y,1)
 
      if (associated(domain2d_spec_out%y)) deallocate(domain2d_spec_out%y) !< Check if allocated
      allocate(domain2d_spec_out%y(starting:ending))
      do i = starting, ending
         call mpp_copy_domain1d_spec(domain2d_spec_in%y(i), domain2d_spec_out%y(i))
      enddo
  endif
 
  end subroutine mpp_copy_domain2d_spec
 
  !> @brief Copies input 1d domain spec to the output 1d domain spec
  subroutine mpp_copy_domain1d_spec(domain1D_spec_in, domain1D_spec_out)
     type(domain1d_spec), intent(in) :: domain1D_spec_in   !< Input
     type(domain1d_spec), intent(out) :: domain1D_spec_out !< Output
 
     domain1d_spec_out%pos = domain1d_spec_in%pos
 
     call mpp_copy_domain_axis_spec(domain1d_spec_in%compute, domain1d_spec_out%compute)
     call mpp_copy_domain_axis_spec(domain1d_spec_in%global, domain1d_spec_out%global)
  end subroutine mpp_copy_domain1d_spec
 
  !> @brief Copies input domain_axis_spec to the output domain_axis_spec
  subroutine mpp_copy_domain_axis_spec(domain_axis_spec_in, domain_axis_spec_out)
     type(domain_axis_spec), intent(in) :: domain_axis_spec_in   !< Input
     type(domain_axis_spec), intent(out) :: domain_axis_spec_out !< Output
 
     domain_axis_spec_out%begin = domain_axis_spec_in%begin
     domain_axis_spec_out%end = domain_axis_spec_in%end
     domain_axis_spec_out%size = domain_axis_spec_in%size
     domain_axis_spec_out%max_size = domain_axis_spec_in%max_size
     domain_axis_spec_out%is_global = domain_axis_spec_in%is_global
  end subroutine mpp_copy_domain_axis_spec
 
  !######################################################################
  subroutine set_group_update(group, domain)
     type(mpp_group_update_type), intent(inout) :: group
     type(domain2d),              intent(inout) :: domain
     integer   :: nscalar, nvector, nlist
     integer   :: nsend, nrecv, nsend_old, nrecv_old
     integer   :: nsend_s, nsend_x, nsend_y
     integer   :: nrecv_s, nrecv_x, nrecv_y
     integer   :: update_buffer_pos, tot_recv_size, tot_send_size
     integer   :: msgsize_s, msgsize_x, msgsize_y, msgsize
     logical   :: recv_s(8), send_s(8)
     logical   :: recv_x(8), send_x(8), recv_y(8), send_y(8)
     integer   :: ntot, n, l, m, ksize
     integer   :: i_s, i_x, i_y, rank_s, rank_x, rank_y, rank
     integer   :: ind_s(3*MAXOVERLAP)
     integer   :: ind_x(3*MAXOVERLAP)
     integer   :: ind_y(3*MAXOVERLAP)
     integer   :: pelist(3*MAXOVERLAP)
     integer   :: send_size(3*MAXOVERLAP)
     integer   :: position_x, position_y, npack, nunpack, dir
     integer   :: pack_buffer_pos, unpack_buffer_pos
     integer   :: omp_get_num_threads, nthreads
     character(len=8)            :: text
     type(overlap_type), pointer :: overPtr => null()
     type(overlapspec), pointer  :: update_s => null()
     type(overlapspec), pointer  :: update_x => null()
     type(overlapspec), pointer  :: update_y => null()
 
     nscalar = group%nscalar
     nvector = group%nvector
 
     !--- get the overlap data type
     select case(group%gridtype)
     case (agrid)
        position_x = center
        position_y = center
     case (bgrid_ne, bgrid_sw)
        position_x = corner
        position_y = corner
     case (cgrid_ne, cgrid_sw)
        position_x = east
        position_y = north
     case (dgrid_ne, dgrid_sw)
        position_x = north
        position_y = east
     case default
        call mpp_error(fatal, "set_group_update: invalid value of gridtype")
     end select
     if(nscalar>0) then
        update_s => search_update_overlap(domain, group%whalo_s, group%ehalo_s, &
             group%shalo_s, group%nhalo_s, group%position)
     endif
     if(nvector>0) then
        update_x => search_update_overlap(domain, group%whalo_v, group%ehalo_v, &
             group%shalo_v, group%nhalo_v, position_x)
        update_y => search_update_overlap(domain, group%whalo_v, group%ehalo_v, &
             group%shalo_v, group%nhalo_v, position_y)
     endif
 
     if(nscalar > 0) then
        recv_s = group%recv_s
        send_s = recv_s
     endif
     if(nvector > 0) then
        recv_x = group%recv_x
        send_x = recv_x
        recv_y = group%recv_y
        send_y = recv_y
     end if
     nlist   = size(domain%list(:))
     group%initialized = .true.
     nsend_s = 0; nsend_x = 0; nsend_y = 0
     nrecv_s = 0; nrecv_x = 0; nrecv_y = 0
 
     if(nscalar > 0) then
         !--- This check could not be done because of memory domain
!        if( group%isize_s .NE. (group%ie_s-group%is_s+1) .OR. group%jsize_s .NE. (group%je_s-group%js_s+1))  &
!           call mpp_error(FATAL, "set_group_update: mismatch of size of the field and domain memory domain")
        nsend_s = update_s%nsend
        nrecv_s = update_s%nrecv
     endif
 
     !--- ksize_s must equal ksize_v
     if(nvector > 0 .AND. nscalar > 0) then
        if(group%ksize_s .NE. group%ksize_v) then
           call mpp_error(fatal, "set_group_update: ksize_s and ksize_v are not equal")
        endif
        ksize = group%ksize_s
     else if (nscalar > 0) then
        ksize = group%ksize_s
     else if (nvector > 0) then
        ksize = group%ksize_v
     else
        call mpp_error(fatal, "set_group_update: nscalar and nvector are all 0")
     endif
 
    nthreads = 1
!$OMP PARALLEL
!$   nthreads = omp_get_num_threads()
!$OMP END PARALLEL
     if( nthreads > nthread_control_loop ) then
        group%k_loop_inside = .false.
     else
        group%k_loop_inside = .true.
     endif
 
     if(nvector > 0) then
        !--- This check could not be done because of memory domain
!        if( group%isize_x .NE. (group%ie_x-group%is_x+1) .OR. group%jsize_x .NE. (group%je_x-group%js_x+1))  &
!           call mpp_error(FATAL, "set_group_update: mismatch of size of the fieldx and domain memory domain")
!        if( group%isize_y .NE. (group%ie_y-group%is_y+1) .OR. group%jsize_y .NE. (group%je_y-group%js_y+1))  &
!           call mpp_error(FATAL, "set_group_update: mismatch of size of the fieldy and domain memory domain")
        nsend_x = update_x%nsend
        nrecv_x = update_x%nrecv
        nsend_y = update_y%nsend
        nrecv_y = update_y%nrecv
     endif
 
     !figure out message size for each processor.
     ntot = nrecv_s + nrecv_x + nrecv_y
     if(ntot > 3*maxoverlap) call mpp_error(fatal, "set_group_update: ntot is greater than 3*MAXOVERLAP")
     n = 1
     i_s = 1
     i_x = 1
     i_y = 1
     ind_s = -1
     ind_x = -1
     ind_y = -1
     nrecv = 0
     do while(n<=ntot)
        if( i_s <= nrecv_s ) then
           rank_s = update_s%recv(i_s)%pe-domain%pe
           if(rank_s .LE. 0) rank_s = rank_s + nlist
        else
           rank_s = -1
        endif
        if( i_x <= nrecv_x ) then
           rank_x = update_x%recv(i_x)%pe-domain%pe
           if(rank_x .LE. 0) rank_x = rank_x + nlist
        else
           rank_x = -1
        endif
        if( i_y <= nrecv_y ) then
           rank_y = update_y%recv(i_y)%pe-domain%pe
           if(rank_y .LE. 0) rank_y = rank_y + nlist
        else
           rank_y = -1
        endif
        nrecv = nrecv + 1
        rank = maxval((/rank_s, rank_x, rank_y/))
        if(rank == rank_s) then
           n = n + 1
           ind_s(nrecv) = i_s
           pelist(nrecv) = update_s%recv(i_s)%pe
           i_s = i_s + 1
        endif
        if(rank == rank_x) then
           n = n + 1
           ind_x(nrecv) = i_x
           pelist(nrecv) = update_x%recv(i_x)%pe
           i_x = i_x + 1
        endif
        if(rank == rank_y) then
           n = n + 1
           ind_y(nrecv) = i_y
           pelist(nrecv) = update_y%recv(i_y)%pe
           i_y = i_y + 1
        endif
     enddo
 
     nrecv_old = nrecv
     nrecv     = 0
     update_buffer_pos = 0
     tot_recv_size = 0
 
     !--- setup for recv
     do l = 1, nrecv_old
        msgsize_s = 0
        msgsize_x = 0
        msgsize_y = 0
        m = ind_s(l)
        if(m>0) msgsize_s = get_mesgsize(update_s%recv(m), recv_s)*ksize*nscalar
        m = ind_x(l)
        if(m>0) msgsize_x = get_mesgsize(update_x%recv(m), recv_x)*ksize*nvector
        m = ind_y(l)
        if(m>0) msgsize_y = get_mesgsize(update_y%recv(m), recv_y)*ksize*nvector
        msgsize = msgsize_s + msgsize_x + msgsize_y
        if( msgsize.GT.0 )then
           tot_recv_size = tot_recv_size + msgsize
           nrecv = nrecv + 1
           if(nrecv > maxoverlap) then
              call mpp_error(fatal, "set_group_update: nrecv is greater than MAXOVERLAP, increase MAXOVERLAP")
           endif
           group%from_pe(nrecv) = pelist(l)
           group%recv_size(nrecv) = msgsize
           group%buffer_pos_recv(nrecv) = update_buffer_pos
           update_buffer_pos = update_buffer_pos + msgsize
        end if
     end do
     group%nrecv = nrecv
 
     !--- setup for unpack
     nunpack = 0
     unpack_buffer_pos = 0
     do l = 1, nrecv_old
        m = ind_s(l)
        if(m>0) then
           overptr => update_s%recv(m)
           do n = 1, overptr%count
              dir = overptr%dir(n)
              if(recv_s(dir)) then
                 nunpack = nunpack + 1
                 if(nunpack > maxoverlap) call mpp_error(fatal,  &
                   "set_group_update: nunpack is greater than MAXOVERLAP, increase MAXOVERLAP 1")
                 group%unpack_type(nunpack) = field_s
                 group%unpack_buffer_pos(nunpack) = unpack_buffer_pos
                 group%unpack_rotation(nunpack) = overptr%rotation(n)
                 group%unpack_is(nunpack) = overptr%is(n)
                 group%unpack_ie(nunpack) = overptr%ie(n)
                 group%unpack_js(nunpack) = overptr%js(n)
                 group%unpack_je(nunpack) = overptr%je(n)
                 group%unpack_size(nunpack) = overptr%msgsize(n)*nscalar
                 unpack_buffer_pos = unpack_buffer_pos + group%unpack_size(nunpack)*ksize
              end if
           end do
        end if
 
        m = ind_x(l)
        if(m>0) then
           overptr => update_x%recv(m)
           do n = 1, overptr%count
              dir = overptr%dir(n)
              if(recv_x(dir)) then
                 nunpack = nunpack + 1
                 if(nunpack > maxoverlap) call mpp_error(fatal,  &
                   "set_group_update: nunpack is greater than MAXOVERLAP, increase MAXOVERLAP 2")
                 group%unpack_type(nunpack) = field_x
                 group%unpack_buffer_pos(nunpack) = unpack_buffer_pos
                 group%unpack_rotation(nunpack) = overptr%rotation(n)
                 group%unpack_is(nunpack) = overptr%is(n)
                 group%unpack_ie(nunpack) = overptr%ie(n)
                 group%unpack_js(nunpack) = overptr%js(n)
                 group%unpack_je(nunpack) = overptr%je(n)
                 group%unpack_size(nunpack) = overptr%msgsize(n)*nvector
                 unpack_buffer_pos = unpack_buffer_pos + group%unpack_size(nunpack)*ksize
              end if
           end do
        end if
 
        m = ind_y(l)
        if(m>0) then
           overptr => update_y%recv(m)
           do n = 1, overptr%count
              dir = overptr%dir(n)
              if(recv_y(dir)) then
                 nunpack = nunpack + 1
                 if(nunpack > maxoverlap) call mpp_error(fatal,  &
                   "set_group_update: nunpack is greater than MAXOVERLAP, increase MAXOVERLAP 3")
                 group%unpack_type(nunpack) = field_y
                 group%unpack_buffer_pos(nunpack) = unpack_buffer_pos
                 group%unpack_rotation(nunpack) = overptr%rotation(n)
                 group%unpack_is(nunpack) = overptr%is(n)
                 group%unpack_ie(nunpack) = overptr%ie(n)
                 group%unpack_js(nunpack) = overptr%js(n)
                 group%unpack_je(nunpack) = overptr%je(n)
                 group%unpack_size(nunpack) = overptr%msgsize(n)*nvector
                 unpack_buffer_pos = unpack_buffer_pos + group%unpack_size(nunpack)*ksize
              end if
           end do
        end if
     end do
     group%nunpack = nunpack
 
     if(update_buffer_pos .NE. unpack_buffer_pos ) call mpp_error(fatal,  &
                   .NE."set_group_update: update_buffer_pos  unpack_buffer_pos")
 
     !figure out message size for each processor.
     ntot = nsend_s + nsend_x + nsend_y
     n = 1
     i_s = 1
     i_x = 1
     i_y = 1
     ind_s = -1
     ind_x = -1
     ind_y = -1
     nsend = 0
     do while(n<=ntot)
        if( i_s <= nsend_s ) then
           rank_s = update_s%send(i_s)%pe-domain%pe
           if(rank_s .LT. 0) rank_s = rank_s + nlist
        else
           rank_s = nlist+1
        endif
        if( i_x <= nsend_x ) then
           rank_x = update_x%send(i_x)%pe-domain%pe
           if(rank_x .LT. 0) rank_x = rank_x + nlist
        else
           rank_x = nlist+1
        endif
        if( i_y <= nsend_y ) then
           rank_y = update_y%send(i_y)%pe-domain%pe
           if(rank_y .LT. 0) rank_y = rank_y + nlist
        else
           rank_y = nlist+1
        endif
        nsend = nsend + 1
        rank = minval((/rank_s, rank_x, rank_y/))
        if(rank == rank_s) then
           n = n + 1
           ind_s(nsend) = i_s
           pelist(nsend) = update_s%send(i_s)%pe
           i_s = i_s + 1
        endif
        if(rank == rank_x) then
           n = n + 1
           ind_x(nsend) = i_x
           pelist(nsend) = update_x%send(i_x)%pe
           i_x = i_x + 1
        endif
        if(rank == rank_y) then
           n = n + 1
           ind_y(nsend) = i_y
           pelist(nsend) = update_y%send(i_y)%pe
           i_y = i_y + 1
        endif
     enddo
 
     nsend_old = nsend
     nsend     = 0
     tot_send_size = 0
     do l = 1, nsend_old
        msgsize_s = 0
        msgsize_x = 0
        msgsize_y = 0
        m = ind_s(l)
        if(m>0) msgsize_s = get_mesgsize(update_s%send(m), send_s)*ksize*nscalar
        m = ind_x(l)
        if(m>0) msgsize_x = get_mesgsize(update_x%send(m), send_x)*ksize*nvector
        m = ind_y(l)
        if(m>0) msgsize_y = get_mesgsize(update_y%send(m), send_y)*ksize*nvector
        msgsize = msgsize_s + msgsize_x + msgsize_y
        if( msgsize.GT.0 )then
           tot_send_size = tot_send_size + msgsize
           nsend = nsend + 1
           if(nsend > maxoverlap) then
              call mpp_error(fatal, "set_group_update: nsend is greater than MAXOVERLAP, increase MAXOVERLAP")
           endif
           send_size(nsend) = msgsize
           group%to_pe(nsend) = pelist(l)
           group%buffer_pos_send(nsend) = update_buffer_pos
           group%send_size(nsend) = msgsize
           update_buffer_pos = update_buffer_pos + msgsize
        end if
     end do
     group%nsend = nsend
 
     !--- setup for pack
     npack = 0
     pack_buffer_pos = unpack_buffer_pos
     do l = 1, nsend_old
        m = ind_s(l)
        if(m>0) then
           overptr => update_s%send(m)
           do n = 1, overptr%count
              dir = overptr%dir(n)
              if(send_s(dir)) then
                 npack = npack + 1
                 if(npack > maxoverlap) call mpp_error(fatal,  &
                   "set_group_update: npack is greater than MAXOVERLAP, increase MAXOVERLAP 1")
                 group%pack_type(npack) = field_s
                 group%pack_buffer_pos(npack) = pack_buffer_pos
                 group%pack_rotation(npack) = overptr%rotation(n)
                 group%pack_is(npack) = overptr%is(n)
                 group%pack_ie(npack) = overptr%ie(n)
                 group%pack_js(npack) = overptr%js(n)
                 group%pack_je(npack) = overptr%je(n)
                 group%pack_size(npack) = overptr%msgsize(n)*nscalar
                 pack_buffer_pos = pack_buffer_pos + group%pack_size(npack)*ksize
              end if
           end do
        end if
 
        m = ind_x(l)
        if(m>0) then
           overptr => update_x%send(m)
           do n = 1, overptr%count
              dir = overptr%dir(n)
              !--- nonsym_edge update is not for rotation of 90 or -90 degree ( cubic sphere grid )
              if( group%nonsym_edge .and. (overptr%rotation(n)==ninety .or. &
                                          overptr%rotation(n)==minus_ninety) ) then
                 call mpp_error(fatal, 'set_group_update: flags=NONSYMEDGEUPDATE is not compatible '// &
                                       'with 90 or -90 degree rotation (normally cubic sphere grid' )
              endif
              if(send_x(dir)) then
                 npack = npack + 1
                 if(npack > maxoverlap) call mpp_error(fatal,  &
                   "set_group_update: npack is greater than MAXOVERLAP, increase MAXOVERLAP 2")
                 group%pack_type(npack) = field_x
                 group%pack_buffer_pos(npack) = pack_buffer_pos
                 group%pack_rotation(npack) = overptr%rotation(n)
                 group%pack_is(npack) = overptr%is(n)
                 group%pack_ie(npack) = overptr%ie(n)
                 group%pack_js(npack) = overptr%js(n)
                 group%pack_je(npack) = overptr%je(n)
                 group%pack_size(npack) = overptr%msgsize(n)*nvector
                 pack_buffer_pos = pack_buffer_pos + group%pack_size(npack)*ksize
              end if
           end do
        end if
 
        m = ind_y(l)
        if(m>0) then
           overptr => update_y%send(m)
           do n = 1, overptr%count
              dir = overptr%dir(n)
              if( group%nonsym_edge .and. (overptr%rotation(n)==ninety .or. &
                                          overptr%rotation(n)==minus_ninety) ) then
                 call mpp_error(fatal, 'set_group_update: flags=NONSYMEDGEUPDATE is not compatible '// &
                                       'with 90 or -90 degree rotation (normally cubic sphere grid' )
              endif
              if(send_y(dir)) then
                 npack = npack + 1
                 if(npack > maxoverlap) call mpp_error(fatal,  &
                   "set_group_update: npack is greater than MAXOVERLAP, increase MAXOVERLAP 3")
                 group%pack_type(npack) = field_y
                 group%pack_buffer_pos(npack) = pack_buffer_pos
                 group%pack_rotation(npack) = overptr%rotation(n)
                 group%pack_is(npack) = overptr%is(n)
                 group%pack_ie(npack) = overptr%ie(n)
                 group%pack_js(npack) = overptr%js(n)
                 group%pack_je(npack) = overptr%je(n)
                 group%pack_size(npack) = overptr%msgsize(n)*nvector
                 pack_buffer_pos = pack_buffer_pos + group%pack_size(npack)*ksize
              end if
           end do
        end if
     end do
     group%npack = npack
     if(update_buffer_pos .NE. pack_buffer_pos ) call mpp_error(fatal,  &
                   .NE."set_group_update: update_buffer_pos  pack_buffer_pos")
 
     !--- make sure the buffer is large enough
     mpp_domains_stack_hwm = max( mpp_domains_stack_hwm, tot_recv_size+tot_send_size )
 
     if( mpp_domains_stack_hwm.GT.mpp_domains_stack_size )then
        write( text,'(i8)' )mpp_domains_stack_hwm
        call mpp_error( fatal, 'set_group_update: mpp_domains_stack overflow, '// &
                     'call mpp_domains_set_stack_size('//trim(text)//') from all PEs.' )
     end if
 
     group%tot_msgsize = tot_recv_size+tot_send_size
 
end subroutine set_group_update
 
 
!######################################################################
  subroutine mpp_clear_group_update(group)
     type(mpp_group_update_type), intent(inout) :: group
 
     group%nscalar = 0
     group%nvector = 0
     group%nsend   = 0
     group%nrecv   = 0
     group%npack   = 0
     group%nunpack = 0
     group%initialized = .false.
 
  end subroutine mpp_clear_group_update
 
!#####################################################################
  function mpp_group_update_initialized(group)
     type(mpp_group_update_type), intent(in) :: group
     logical :: mpp_group_update_initialized
 
     mpp_group_update_initialized = group%initialized
 
  end function mpp_group_update_initialized
 
!#####################################################################
  function mpp_group_update_is_set(group)
     type(mpp_group_update_type), intent(in) :: group
     logical :: mpp_group_update_is_set
 
     mpp_group_update_is_set = (group%nscalar > 0 .OR. group%nvector > 0)
 
  end function mpp_group_update_is_set
!> @}