mosaic/grid.F90

   1 !***********************************************************************
   2 !*                   GNU Lesser General Public License
   3 !*
   4 !* This file is part of the GFDL Flexible Modeling System (FMS).
   5 !*
   6 !* FMS is free software: you can redistribute it and/or modify it under
   7 !* the terms of the GNU Lesser General Public License as published by
   8 !* the Free Software Foundation, either version 3 of the License, or (at
   9 !* your option) any later version.
  10 !*
  11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
  12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 !* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 !* for more details.
  15 !*
  16 !* You should have received a copy of the GNU Lesser General Public
  17 !* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
  18 !***********************************************************************
  19 !> @defgroup grid_mod grid_mod
  20 !> @ingroup mosaic
  21 !> @brief Routines for grid calculations
  22
  23 module grid_mod
  24 #ifdef use_deprecated_io
  25
  26 use mpp_mod, only : mpp_root_pe, uppercase, lowercase, FATAL, NOTE, mpp_error
  27 use constants_mod, only : PI, radius
  28 use fms_io_mod, only : get_great_circle_algorithm, get_global_att_value, string, &
  29                        field_exist, field_size, read_data
  30 use mosaic_mod, only : get_mosaic_ntiles, get_mosaic_xgrid_size, get_mosaic_grid_sizes, &
  31      get_mosaic_xgrid, calc_mosaic_grid_area, calc_mosaic_grid_great_circle_area
  32
  33 ! the following two use statement are only needed for define_cube_mosaic
  34 use mpp_domains_mod, only : domain2d, mpp_define_mosaic, mpp_get_compute_domain, &
  35                             mpp_get_global_domain, domainUG, mpp_pass_SG_to_UG
  36 use mosaic_mod, only : get_mosaic_ncontacts, get_mosaic_contact
  37
  38 implicit none;private
  39
  40 ! ==== public interfaces =====================================================
  41 ! grid dimension inquiry subroutines
  42 public :: get_grid_ntiles
  43 public :: get_grid_size
  44 ! grid geometry inquiry subroutines
  45 public :: get_grid_cell_centers
  46 public :: get_grid_cell_vertices
  47 ! grid area inquiry subroutines
  48 public :: get_grid_cell_area
  49 public :: get_grid_comp_area
  50 ! decompose cubed sphere domains -- probably does not belong here, but it should
  51 ! be in some place available for component models
  52 public :: define_cube_mosaic
  53 ! ==== end of public interfaces ==============================================
  54
  55 !> returns horizontal sizes of the grid
  56 !> @ingroup grid_mod
  57 interface get_grid_size
  58    module procedure get_grid_size_for_all_tiles
  59    module procedure get_grid_size_for_one_tile
  60 end interface
  61 !> returns number of tiles
  62 !> @ingroup grid_mod
  63 interface get_grid_cell_vertices
  64    module procedure get_grid_cell_vertices_1D
  65    module procedure get_grid_cell_vertices_2D
  66    module procedure get_grid_cell_vertices_UG
  67 end interface
  68
  69 !> @ingroup grid_mod
  70 interface get_grid_cell_centers
  71    module procedure get_grid_cell_centers_1D
  72    module procedure get_grid_cell_centers_2D
  73    module procedure get_grid_cell_centers_UG
  74 end interface
  75
  76 !> @ingroup grid_mod
  77 interface get_grid_cell_area
  78    module procedure get_grid_cell_area_SG
  79    module procedure get_grid_cell_area_UG
  80 end interface get_grid_cell_area
  81
  82 !> @ingroup grid_mod
  83 interface get_grid_comp_area
  84    module procedure get_grid_comp_area_SG
  85    module procedure get_grid_comp_area_UG
  86 end interface get_grid_comp_area
  87
  88 !> @addtogroup grid_mod
  89 !> @{
  90
  91 ! ==== module constants ======================================================
  92 character(len=*), parameter :: &
  93      module_name = 'grid_mod'
  94
  95 ! Include variable "version" to be written to log file.
  96 #include<file_version.h>
  97
  98 character(len=*), parameter :: &
  99      grid_dir  = 'INPUT/',     &      !< root directory for all grid files
 100      grid_file = 'INPUT/grid_spec.nc' !< name of the grid spec file
 101
 102 integer, parameter :: &
 103      MAX_NAME = 256,  & !< max length of the variable names
 104      MAX_FILE = 1024, & !< max length of the file names
 105      VERSION_0 = 0,   &
 106      VERSION_1 = 1,   &
 107      VERSION_2 = 2
 108
 109 integer, parameter :: BUFSIZE = 1048576  !< This is used to control memory usage in get_grid_comp_area
 110                                          !! We may change this to a namelist variable is needed.
 111
 112 ! ==== module variables ======================================================
 113 integer :: grid_version = -1
 114 logical :: great_circle_algorithm = .FALSE.
 115 logical :: first_call = .TRUE.
 116
 117
 118 contains
 119
 120 function get_grid_version()
 121   integer :: get_grid_version
 122
 123   if(first_call) then
 124      great_circle_algorithm = get_great_circle_algorithm()
 125      first_call = .FALSE.
 126   endif
 127
 128   if(grid_version<0) then
 129     if(field_exist(grid_file, 'geolon_t')) then
 130        grid_version = VERSION_0
 131     else if(field_exist(grid_file, 'x_T')) then
 132        grid_version = VERSION_1
 133     else if(field_exist(grid_file, 'ocn_mosaic_file') ) then
 134        grid_version = VERSION_2
 135     else
 136        call mpp_error(FATAL, module_name//&
 137                     & '/get_grid_version: Can''t determine the version of the grid spec:'// &
 138                     & ' none of "x_T", "geolon_t", or "ocn_mosaic_file" exist in file "'//trim(grid_file)//'"')
 139     endif
 140   endif
 141   get_grid_version = grid_version
 142 end function get_grid_version
 143
 144
 145 ! ============================================================================
 146 ! ============================================================================
 147 !> Returns number of tiles for a given component
 148 subroutine get_grid_ntiles(component,ntiles)
 149   character(len=*)     :: component
 150   integer, intent(out) :: ntiles
 151
 152   ! local vars
 153   character(len=MAX_FILE) :: component_mosaic
 154
 155   select case (get_grid_version())
 156   case(VERSION_0,VERSION_1)
 157      ntiles = 1
 158   case(VERSION_2)
 159      call read_data(grid_file,trim(lowercase(component))//'_mosaic_file',component_mosaic)
 160      ntiles = get_mosaic_ntiles(grid_dir//trim(component_mosaic))
 161   end select
 162 end subroutine get_grid_ntiles
 163
 164
 165 ! ============================================================================
 166 ! ============================================================================
 167 !> Returns size of the grid for each of the tiles
 168 subroutine get_grid_size_for_all_tiles(component,nx,ny)
 169   character(len=*)     :: component
 170   integer, intent(inout) :: nx(:),ny(:)
 171
 172   ! local vars
 173   integer :: siz(4) ! for the size of external fields
 174   character(len=MAX_NAME) :: varname1, varname2
 175   character(len=MAX_FILE) :: component_mosaic
 176
 177   varname1 = 'AREA_'//trim(uppercase(component))
 178   varname2 = trim(lowercase(component))//'_mosaic_file'
 179
 180   select case (get_grid_version())
 181   case(VERSION_0,VERSION_1)
 182      call field_size(grid_file, varname1, siz)
 183      nx(1) = siz(1); ny(1)=siz(2)
 184   case(VERSION_2) ! mosaic file
 185      call read_data(grid_file,varname2, component_mosaic)
 186      call get_mosaic_grid_sizes(grid_dir//trim(component_mosaic),nx,ny)
 187   end select
 188 end subroutine get_grid_size_for_all_tiles
 189
 190
 191 ! ============================================================================
 192 ! ============================================================================
 193 !> Returns size of the grid for one of the tiles
 194 subroutine get_grid_size_for_one_tile(component,tile,nx,ny)
 195   character(len=*)       :: component
 196   integer, intent(in)    :: tile
 197   integer, intent(inout) :: nx,ny
 198
 199   ! local vars
 200   integer, allocatable :: nnx(:), nny(:)
 201   integer :: ntiles
 202
 203   call get_grid_ntiles(component, ntiles)
 204   if(tile>0.and.tile<=ntiles) then
 205      allocate(nnx(ntiles),nny(ntiles))
 206      call get_grid_size_for_all_tiles(component,nnx,nny)
 207      nx = nnx(tile); ny = nny(tile)
 208      deallocate(nnx,nny)
 209   else
 210      call mpp_error(FATAL, 'get_grid_size: requested tile index '// &
 211                     & trim(string(tile))//' is out of bounds (1:'//trim(string(ntiles))//')')
 212   endif
 213 end subroutine get_grid_size_for_one_tile
 214
 215 ! ============================================================================
 216 ! ============================================================================
 217 !> Return grid cell area for the specified model component and tile
 218 subroutine get_grid_cell_area_SG(component, tile, cellarea, domain)
 219   character(len=*), intent(in)    :: component
 220   integer         , intent(in)    :: tile
 221   real            , intent(inout) :: cellarea(:,:)
 222   type(domain2d)  , intent(in), optional :: domain
 223
 224   ! local vars
 225   integer :: nlon, nlat
 226   real, allocatable :: glonb(:,:), glatb(:,:)
 227
 228   select case(get_grid_version())
 229   case(VERSION_0,VERSION_1)
 230      select case(trim(component))
 231      case('LND')
 232         call read_data(grid_file, 'AREA_LND_CELL', cellarea, &
 233             no_domain=.not.present(domain), domain=domain)
 234      case('ATM','OCN')
 235         call read_data(grid_file, 'AREA_'//trim(uppercase(component)),cellarea,&
 236             no_domain=.not.present(domain),domain=domain)
 237      case default
 238         call mpp_error(FATAL, module_name//'/get_grid_cell_area: Illegal component name "'//trim(component) &
 239                             & //'": must be one of ATM, LND, or OCN')
 240      end select
 241      ! convert area to m2
 242      cellarea = cellarea*4.*PI*radius**2
 243   case(VERSION_2)
 244      if (present(domain)) then
 245         call mpp_get_compute_domain(domain,xsize=nlon,ysize=nlat)
 246      else
 247         call get_grid_size(component,tile,nlon,nlat)
 248      endif
 249      allocate(glonb(nlon+1,nlat+1),glatb(nlon+1,nlat+1))
 250      call get_grid_cell_vertices(component, tile, glonb, glatb, domain)
 251      if (great_circle_algorithm) then
 252         call calc_mosaic_grid_great_circle_area(glonb*pi/180.0, glatb*pi/180.0, cellarea)
 253      else
 254         call calc_mosaic_grid_area(glonb*pi/180.0, glatb*pi/180.0, cellarea)
 255      end if
 256      deallocate(glonb,glatb)
 257   end select
 258
 259 end subroutine get_grid_cell_area_SG
 260
 261 ! ============================================================================
 262 ! ============================================================================
 263 !> Get the area of the component per grid cell
 264 subroutine get_grid_comp_area_SG(component,tile,area,domain)
 265   character(len=*) :: component
 266   integer, intent(in) :: tile
 267   real, intent(inout) :: area(:,:)
 268   type(domain2d), intent(in), optional :: domain
 269   ! local vars
 270   integer :: n_xgrid_files ! number of exchange grid files in the mosaic
 271   integer :: siz(4), nxgrid
 272   integer :: i,j,m,n
 273   integer, allocatable :: i1(:), j1(:), i2(:), j2(:)
 274   real, allocatable :: xgrid_area(:)
 275   real, allocatable :: rmask(:,:)
 276   character(len=MAX_NAME) :: &
 277      xgrid_name, & ! name of the variable holding xgrid names
 278      tile_name,  & ! name of the tile
 279      xgrid_file, & ! name of the current xgrid file
 280      mosaic_name,& ! name of the mosaic
 281      mosaic_file,&
 282      tilefile
 283   character(len=4096)     :: attvalue
 284   character(len=MAX_NAME), allocatable :: nest_tile_name(:)
 285   integer :: is,ie,js,je ! boundaries of our domain
 286   integer :: i0, j0 ! offsets for x and y, respectively
 287   integer :: num_nest_tile, ntiles
 288   logical :: is_nest
 289   integer :: found_xgrid_files ! how many xgrid files we actually found in the grid spec
 290   integer :: ibegin, iend, bsize, l
 291
 292   select case (get_grid_version())
 293   case(VERSION_0,VERSION_1)
 294      select case(component)
 295      case('ATM')
 296         call read_data(grid_file,'AREA_ATM',area, no_domain=.not.present(domain),domain=domain)
 297      case('OCN')
 298         allocate(rmask(size(area,1),size(area,2)))
 299         call read_data(grid_file,'AREA_OCN',area, no_domain=.not.present(domain),domain=domain)
 300         call read_data(grid_file,'wet',     rmask,no_domain=.not.present(domain),domain=domain)
 301         area = area*rmask
 302         deallocate(rmask)
 303      case('LND')
 304         call read_data(grid_file,'AREA_LND',area,no_domain=.not.present(domain),domain=domain)
 305      case default
 306         call mpp_error(FATAL, module_name// &
 307               & '/get_grid_comp_area: Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN')
 308      end select
 309   case(VERSION_2) ! mosaic gridspec
 310      select case (component)
 311      case ('ATM')
 312         ! just read the grid cell area and return
 313         call get_grid_cell_area(component,tile,area)
 314         return
 315      case ('LND')
 316         xgrid_name = 'aXl_file'
 317         call read_data(grid_file, 'lnd_mosaic', mosaic_name)
 318         tile_name  = trim(mosaic_name)//'_tile'//char(tile+ichar('0'))
 319      case ('OCN')
 320         xgrid_name = 'aXo_file'
 321         call read_data(grid_file, 'ocn_mosaic', mosaic_name)
 322         tile_name  = trim(mosaic_name)//'_tile'//char(tile+ichar('0'))
 323      case default
 324         call mpp_error(FATAL, module_name// &
 325               & '/get_grid_comp_area: Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN')
 326      end select
 327      ! get the boundaries of the requested domain
 328      if(present(domain)) then
 329         call mpp_get_compute_domain(domain,is,ie,js,je)
 330         i0 = 1-is ; j0=1-js
 331      else
 332         call get_grid_size(component,tile,ie,je)
 333         is = 1 ; i0 = 0
 334         js = 1 ; j0 = 0
 335      endif
 336      if (size(area,1)/=ie-is+1.or.size(area,2)/=je-js+1) &
 337         call mpp_error(FATAL, module_name// &
 338                        & '/get_grid_comp_area: size of the output argument "area" is not consistent with the domain')
 339
 340      ! find the nest tile
 341      call read_data(grid_file, 'atm_mosaic', mosaic_name)
 342      call read_data(grid_file,'atm_mosaic_file',mosaic_file)
 343      mosaic_file = grid_dir//trim(mosaic_file)
 344      ntiles = get_mosaic_ntiles(trim(mosaic_file))
 345      allocate(nest_tile_name(ntiles))
 346      num_nest_tile = 0
 347      do n = 1, ntiles
 348         call read_data(mosaic_file, 'gridfiles', tilefile, level=n)
 349         tilefile = grid_dir//trim(tilefile)
 350         if( get_global_att_value(tilefile, "nest_grid", attvalue) ) then
 351            if(trim(attvalue) == "TRUE") then
 352               num_nest_tile = num_nest_tile + 1
 353               nest_tile_name(num_nest_tile) = trim(mosaic_name)//'_tile'//char(n+ichar('0'))
 354            else if(trim(attvalue) .NE. "FALSE") then
 355               call mpp_error(FATAL, module_name//'/get_grid_comp_area: value of global attribute nest_grid in file'// &
 356                              &  trim(tilefile)//' should be TRUE of FALSE')
 357            endif
 358         end if
 359      end do
 360      area(:,:) = 0.
 361      if(field_exist(grid_file,xgrid_name)) then
 362         ! get the number of the exchange-grid files
 363         call field_size(grid_file,xgrid_name,siz)
 364         n_xgrid_files = siz(2)
 365         found_xgrid_files = 0
 366         ! loop through all exchange grid files
 367         do n = 1, n_xgrid_files
 368            ! get the name of the current exchange grid file
 369            call read_data(grid_file,xgrid_name,xgrid_file,level=n)
 370            ! skip the rest of the loop if the name of the current tile isn't found
 371            ! in the file name, but check this only if there is more than 1 tile
 372            if(n_xgrid_files>1) then
 373               if(index(xgrid_file,trim(tile_name))==0) cycle
 374            endif
 375            found_xgrid_files = found_xgrid_files + 1
 376            !---make sure the atmosphere grid is not a nested grid
 377            is_nest = .false.
 378            do m = 1, num_nest_tile
 379               if(index(xgrid_file, trim(nest_tile_name(m))) .NE. 0) then
 380                  is_nest = .true.
 381                  exit
 382               end if
 383            end do
 384            if(is_nest) cycle
 385
 386            ! finally read the exchange grid
 387            nxgrid = get_mosaic_xgrid_size(grid_dir//xgrid_file)
 388            if(nxgrid < BUFSIZE) then
 389               allocate(i1(nxgrid), j1(nxgrid), i2(nxgrid), j2(nxgrid), xgrid_area(nxgrid))
 390            else
 391               allocate(i1(BUFSIZE), j1(BUFSIZE), i2(BUFSIZE), j2(BUFSIZE), xgrid_area(BUFSIZE))
 392            endif
 393            ibegin = 1
 394            do l = 1,nxgrid,BUFSIZE
 395               bsize = min(BUFSIZE, nxgrid-l+1)
 396               iend = ibegin + bsize - 1
 397               call get_mosaic_xgrid(grid_dir//xgrid_file, i1(1:bsize), j1(1:bsize), i2(1:bsize), j2(1:bsize), &
 398                                     xgrid_area(1:bsize), ibegin, iend)
 399               ! and sum the exchange grid areas
 400               do m = 1, bsize
 401                  i = i2(m); j = j2(m)
 402                  if (i<is.or.i>ie) cycle
 403                  if (j<js.or.j>je) cycle
 404                  area(i+i0,j+j0) = area(i+i0,j+j0) + xgrid_area(m)
 405               end do
 406               ibegin = iend + 1
 407            enddo
 408            deallocate(i1, j1, i2, j2, xgrid_area)
 409         enddo
 410         if (found_xgrid_files == 0) &
 411            call mpp_error(FATAL, 'get_grid_comp_area: no xgrid files were found for component '// &
 412                           & trim(component)//' (mosaic name is '//trim(mosaic_name)//')')
 413
 414      endif
 415      deallocate(nest_tile_name)
 416   end select ! version
 417   ! convert area to m2
 418   area = area*4.*PI*radius**2
 419 end subroutine get_grid_comp_area_SG
 420
 421 !======================================================================
 422 subroutine get_grid_cell_area_UG(component, tile, cellarea, SG_domain, UG_domain)
 423   character(len=*),   intent(in)    :: component
 424   integer         ,   intent(in)    :: tile
 425   real            ,   intent(inout) :: cellarea(:)
 426   type(domain2d)  ,   intent(in)    :: SG_domain
 427   type(domainUG)  ,   intent(in)    :: UG_domain
 428   integer :: is, ie, js, je
 429   real, allocatable :: SG_area(:,:)
 430
 431   call mpp_get_compute_domain(SG_domain, is, ie, js, je)
 432   allocate(SG_area(is:ie, js:je))
 433   call get_grid_cell_area_SG(component, tile, SG_area, SG_domain)
 434   call mpp_pass_SG_to_UG(UG_domain, SG_area, cellarea)
 435   deallocate(SG_area)
 436
 437 end subroutine get_grid_cell_area_UG
 438
 439 subroutine get_grid_comp_area_UG(component, tile, area, SG_domain, UG_domain)
 440   character(len=*),   intent(in)    :: component
 441   integer         ,   intent(in)    :: tile
 442   real            ,   intent(inout) :: area(:)
 443   type(domain2d)  ,   intent(in)    :: SG_domain
 444   type(domainUG)  ,   intent(in)    :: UG_domain
 445   integer :: is, ie, js, je
 446   real, allocatable :: SG_area(:,:)
 447
 448   call mpp_get_compute_domain(SG_domain, is, ie, js, je)
 449   allocate(SG_area(is:ie, js:je))
 450   call get_grid_comp_area_SG(component, tile, SG_area, SG_domain)
 451   call mpp_pass_SG_to_UG(UG_domain, SG_area, area)
 452   deallocate(SG_area)
 453
 454 end subroutine get_grid_comp_area_UG
 455
 456
 457 ! ============================================================================
 458 ! ============================================================================
 459 !> Returns arrays of global grid cell boundaries for given model component and
 460 !! mosaic tile number.
 461 !!
 462 !> @note In the case of non-lat-lon grid the returned coordinates may have be not so
 463 !! meaningful, by the very nature of such grids. But presumably these 1D coordinate
 464 !! arrays are good enough for diag axis and such.
 465 subroutine get_grid_cell_vertices_1D(component, tile, glonb, glatb)
 466   character(len=*), intent(in) :: component
 467   integer,          intent(in) :: tile
 468   real,          intent(inout) :: glonb(:),glatb(:)
 469
 470   integer                      :: nlon, nlat
 471   integer                      :: start(4), nread(4)
 472   real, allocatable            :: tmp(:,:), x_vert_t(:,:,:), y_vert_t(:,:,:)
 473   character(len=MAX_FILE)      :: filename1, filename2
 474
 475   call get_grid_size_for_one_tile(component, tile, nlon, nlat)
 476   if (size(glonb(:))/=nlon+1) &
 477        call mpp_error (FATAL,  module_name// &
 478                        & '/get_grid_cell_vertices_1D: Size of argument "glonb" is not consistent with the grid size')
 479   if (size(glatb(:))/=nlat+1) &
 480        call mpp_error (FATAL, module_name// &
 481                        & '/get_grid_cell_vertices_1D: Size of argument "glatb" is not consistent with the grid size')
 482   if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
 483      call mpp_error(FATAL, module_name//'/get_grid_cell_vertices_1D: Illegal component name "'// &
 484                     & trim(component)//'": must be one of ATM, LND, or OCN')
 485   endif
 486
 487   select case(get_grid_version())
 488   case(VERSION_0)
 489      select case(trim(component))
 490      case('ATM','LND')
 491         call read_data(grid_file, 'xb'//lowercase(component(1:1)), glonb, no_domain=.true.)
 492         call read_data(grid_file, 'yb'//lowercase(component(1:1)), glatb, no_domain=.true.)
 493      case('OCN')
 494         call read_data(grid_file, "gridlon_vert_t", glonb, no_domain=.true.)
 495         call read_data(grid_file, "gridlat_vert_t", glatb, no_domain=.true.)
 496      end select
 497   case(VERSION_1)
 498      select case(trim(component))
 499      case('ATM','LND')
 500         call read_data(grid_file, 'xb'//lowercase(component(1:1)), glonb, no_domain=.true.)
 501         call read_data(grid_file, 'yb'//lowercase(component(1:1)), glatb, no_domain=.true.)
 502      case('OCN')
 503         allocate (x_vert_t(nlon,1,2), y_vert_t(1,nlat,2) )
 504         start = 1; nread = 1
 505         nread(1) = nlon; nread(2) = 1; start(3) = 1
 506         call read_data(grid_file, "x_vert_T", x_vert_t(:,:,1), start, nread, no_domain=.TRUE.)
 507         nread(1) = nlon; nread(2) = 1; start(3) = 2
 508         call read_data(grid_file, "x_vert_T", x_vert_t(:,:,2), start, nread, no_domain=.TRUE.)
 509
 510         nread(1) = 1; nread(2) = nlat; start(3) = 1
 511         call read_data(grid_file, "y_vert_T", y_vert_t(:,:,1), start, nread, no_domain=.TRUE.)
 512         nread(1) = 1; nread(2) = nlat; start(3) = 4
 513         call read_data(grid_file, "y_vert_T", y_vert_t(:,:,2), start, nread, no_domain=.TRUE.)
 514         glonb(1:nlon) = x_vert_t(1:nlon,1,1)
 515         glonb(nlon+1) = x_vert_t(nlon,1,2)
 516         glatb(1:nlat) = y_vert_t(1,1:nlat,1)
 517         glatb(nlat+1) = y_vert_t(1,nlat,2)
 518         deallocate(x_vert_t, y_vert_t)
 519      end select
 520   case(VERSION_2)
 521      ! get the name of the mosaic file for the component
 522      call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
 523      filename1=grid_dir//trim(filename1)
 524      ! get the name of the grid file for the component and tile
 525      call read_data(filename1, 'gridfiles', filename2, level=tile)
 526      filename2 = grid_dir//trim(filename2)
 527
 528      start = 1; nread = 1
 529      nread(1) = 2*nlon+1
 530      allocate( tmp(2*nlon+1,1) )
 531      call read_data(filename2, "x", tmp, start, nread, no_domain=.TRUE.)
 532      glonb(1:nlon+1) = tmp(1:2*nlon+1:2,1)
 533      deallocate(tmp)
 534      allocate(tmp(1,2*nlat+1))
 535
 536      start = 1; nread = 1
 537      nread(2) = 2*nlat+1
 538      call read_data(filename2, "y", tmp, start, nread, no_domain=.TRUE.)
 539      glatb(1:nlat+1) = tmp(1,1:2*nlat+1:2)
 540      deallocate(tmp)
 541   end select
 542
 543 end subroutine get_grid_cell_vertices_1D
 544
 545 ! ============================================================================
 546 ! ============================================================================
 547 !> Returns cell vertices for the specified model component and mosaic tile number
 548 subroutine get_grid_cell_vertices_2D(component, tile, lonb, latb, domain)
 549   character(len=*),         intent(in) :: component
 550   integer,                  intent(in) :: tile
 551   real,                  intent(inout) :: lonb(:,:),latb(:,:)
 552   type(domain2d), optional, intent(in) :: domain
 553
 554   ! local vars
 555   character(len=MAX_FILE) :: filename1, filename2
 556   integer :: nlon, nlat
 557   integer :: i,j
 558   real, allocatable :: buffer(:), tmp(:,:), x_vert_t(:,:,:), y_vert_t(:,:,:)
 559   integer :: is,ie,js,je ! boundaries of our domain
 560   integer :: i0,j0 ! offsets for coordinates
 561   integer :: isg, jsg
 562   integer :: start(4), nread(4)
 563
 564   call get_grid_size_for_one_tile(component, tile, nlon, nlat)
 565   if (present(domain)) then
 566     call mpp_get_compute_domain(domain,is,ie,js,je)
 567   else
 568     is = 1 ; ie = nlon
 569     js = 1 ; je = nlat
 570     !--- domain normally should be present
 571     call mpp_error (NOTE, module_name//'/get_grid_cell_vertices: domain is not present, global data will be read')
 572   endif
 573   i0 = -is+1; j0 = -js+1
 574
 575   ! verify that lonb and latb sizes are consistent with the size of domain
 576   if (size(lonb,1)/=ie-is+2.or.size(lonb,2)/=je-js+2) &
 577        call mpp_error (FATAL, module_name// &
 578                        & '/get_grid_cell_vertices: Size of argument "lonb" is not consistent with the domain size')
 579   if (size(latb,1)/=ie-is+2.or.size(latb,2)/=je-js+2) &
 580        call mpp_error (FATAL, module_name// &
 581                        & '/get_grid_cell_vertices: Size of argument "latb" is not consistent with the domain size')
 582   if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
 583      call mpp_error(FATAL, module_name//'/get_grid_cell_vertices: Illegal component name "'// &
 584                     & trim(component)//'": must be one of ATM, LND, or OCN')
 585   endif
 586
 587   select case(get_grid_version())
 588   case(VERSION_0)
 589      select case(component)
 590      case('ATM','LND')
 591         allocate(buffer(max(nlon,nlat)+1))
 592         ! read coordinates of grid cell vertices
 593         call read_data(grid_file, 'xb'//lowercase(component(1:1)), buffer(1:nlon+1), no_domain=.true.)
 594         do j = js, je+1
 595            do i = is, ie+1
 596               lonb(i+i0,j+j0) = buffer(i)
 597            enddo
 598         enddo
 599         call read_data(grid_file, 'yb'//lowercase(component(1:1)), buffer(1:nlat+1), no_domain=.true.)
 600         do j = js, je+1
 601            do i = is, ie+1
 602               latb(i+i0,j+j0) = buffer(j)
 603            enddo
 604         enddo
 605         deallocate(buffer)
 606      case('OCN')
 607         if (present(domain)) then
 608            start = 1; nread = 1
 609            start(1) = is; start(2) = js
 610            nread(1) = ie-is+2; nread(2) = je-js+2
 611            call read_data(grid_file, 'geolon_vert_t', lonb, start, nread, no_domain=.true. )
 612            call read_data(grid_file, 'geolat_vert_t', latb, start, nread, no_domain=.true. )
 613          else
 614            call read_data(grid_file, 'geolon_vert_t', lonb, no_domain=.TRUE. )
 615            call read_data(grid_file, 'geolat_vert_t', latb, no_domain=.TRUE. )
 616          endif
 617      end select
 618   case(VERSION_1)
 619      select case(component)
 620      case('ATM','LND')
 621         allocate(buffer(max(nlon,nlat)+1))
 622         ! read coordinates of grid cell vertices
 623         call read_data(grid_file, 'xb'//lowercase(component(1:1)), buffer(1:nlon+1), no_domain=.true.)
 624         do j = js, je+1
 625            do i = is, ie+1
 626               lonb(i+i0,j+j0) = buffer(i)
 627            enddo
 628         enddo
 629         call read_data(grid_file, 'yb'//lowercase(component(1:1)), buffer(1:nlat+1), no_domain=.true.)
 630         do j = js, je+1
 631            do i = is, ie+1
 632               latb(i+i0,j+j0) = buffer(j)
 633            enddo
 634         enddo
 635         deallocate(buffer)
 636      case('OCN')
 637         nlon=ie-is+1; nlat=je-js+1
 638         allocate (x_vert_t(nlon,nlat,4), y_vert_t(nlon,nlat,4) )
 639         call read_data(grid_file, 'x_vert_T', x_vert_t, no_domain=.not.present(domain), domain=domain )
 640         call read_data(grid_file, 'y_vert_T', y_vert_t, no_domain=.not.present(domain), domain=domain )
 641         lonb(1:nlon,1:nlat) = x_vert_t(1:nlon,1:nlat,1)
 642         lonb(nlon+1,1:nlat) = x_vert_t(nlon,1:nlat,2)
 643         lonb(1:nlon,nlat+1) = x_vert_t(1:nlon,nlat,4)
 644         lonb(nlon+1,nlat+1) = x_vert_t(nlon,nlat,3)
 645         latb(1:nlon,1:nlat) = y_vert_t(1:nlon,1:nlat,1)
 646         latb(nlon+1,1:nlat) = y_vert_t(nlon,1:nlat,2)
 647         latb(1:nlon,nlat+1) = y_vert_t(1:nlon,nlat,4)
 648         latb(nlon+1,nlat+1) = y_vert_t(nlon,nlat,3)
 649         deallocate(x_vert_t, y_vert_t)
 650      end select
 651   case(VERSION_2)
 652      ! get the name of the mosaic file for the component
 653      call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
 654      filename1=grid_dir//trim(filename1)
 655      ! get the name of the grid file for the component and tile
 656      call read_data(filename1, 'gridfiles', filename2, level=tile)
 657      filename2 = grid_dir//trim(filename2)
 658      if(PRESENT(domain)) then
 659         call mpp_get_global_domain(domain, xbegin=isg, ybegin=jsg)
 660         start = 1; nread = 1
 661         start(1) = 2*(is-isg+1) - 1; nread(1) = 2*(ie-is)+3
 662         start(2) = 2*(js-jsg+1) - 1; nread(2) = 2*(je-js)+3
 663         allocate(tmp(nread(1), nread(2)) )
 664         call read_data(filename2, 'x', tmp, start, nread, no_domain=.TRUE.)
 665         do j = 1, je-js+2
 666            do i = 1, ie-is+2
 667               lonb(i,j) = tmp(2*i-1,2*j-1)
 668            enddo
 669         enddo
 670         call read_data(filename2, 'y', tmp, start, nread, no_domain=.TRUE.)
 671         do j = 1, je-js+2
 672            do i = 1, ie-is+2
 673               latb(i,j) = tmp(2*i-1,2*j-1)
 674            enddo
 675         enddo
 676      else
 677         allocate(tmp(2*nlon+1,2*nlat+1))
 678         call read_data(filename2, 'x', tmp, no_domain=.TRUE.)
 679         do j = js, je+1
 680            do i = is, ie+1
 681               lonb(i+i0,j+j0) = tmp(2*i-1,2*j-1)
 682            end do
 683         end do
 684         call read_data(filename2, 'y', tmp, no_domain=.TRUE.)
 685         do j = js, je+1
 686            do i = is, ie+1
 687               latb(i+i0,j+j0) = tmp(2*i-1,2*j-1)
 688            end do
 689         end do
 690      endif
 691      deallocate(tmp)
 692   end select
 693
 694 end subroutine get_grid_cell_vertices_2D
 695
 696
 697 subroutine get_grid_cell_vertices_UG(component, tile, lonb, latb, SG_domain, UG_domain)
 698   character(len=*),         intent(in) :: component
 699   integer,                  intent(in) :: tile
 700   real,                  intent(inout) :: lonb(:,:),latb(:,:) ! The second dimension is 4
 701   type(domain2d)  ,   intent(in)       :: SG_domain
 702   type(domainUG)  ,   intent(in)       :: UG_domain
 703   integer :: is, ie, js, je, i, j
 704   real, allocatable :: SG_lonb(:,:), SG_latb(:,:), tmp(:,:,:)
 705
 706   call mpp_get_compute_domain(SG_domain, is, ie, js, je)
 707   allocate(SG_lonb(is:ie+1, js:je+1))
 708   allocate(SG_latb(is:ie+1, js:je+1))
 709   allocate(tmp(is:ie,js:je,4))
 710   call get_grid_cell_vertices_2D(component, tile, SG_lonb, SG_latb, SG_domain)
 711   do j = js, je
 712      do i = is, ie
 713         tmp(i,j,1) = SG_lonb(i,j)
 714         tmp(i,j,2) = SG_lonb(i+1,j)
 715         tmp(i,j,3) = SG_lonb(i+1,j+1)
 716         tmp(i,j,4) = SG_lonb(i,j+1)
 717      enddo
 718   enddo
 719   call mpp_pass_SG_to_UG(UG_domain, tmp, lonb)
 720   do j = js, je
 721      do i = is, ie
 722         tmp(i,j,1) = SG_latb(i,j)
 723         tmp(i,j,2) = SG_latb(i+1,j)
 724         tmp(i,j,3) = SG_latb(i+1,j+1)
 725         tmp(i,j,4) = SG_latb(i,j+1)
 726      enddo
 727   enddo
 728   call mpp_pass_SG_to_UG(UG_domain, tmp, latb)
 729
 730
 731   deallocate(SG_lonb, SG_latb, tmp)
 732
 733 end subroutine get_grid_cell_vertices_UG
 734
 735 ! ============================================================================
 736 !> Returns global coordinate arrays fro given model component and mosaic tile number
 737 !! @note In the case of non-lat-lon grid those coordinates may have be not so
 738 !! meaningful, by the very nature of such grids. But presumably these 1D coordinate
 739 !! arrays are good enough for diag axis and such.
 740 subroutine get_grid_cell_centers_1D(component, tile, glon, glat)
 741   character(len=*), intent(in) :: component
 742   integer, intent(in) :: tile
 743   real, intent(inout) :: glon(:),glat(:)
 744   integer                      :: nlon, nlat
 745   integer                      :: start(4), nread(4)
 746   real, allocatable            :: tmp(:,:)
 747   character(len=MAX_FILE)      :: filename1, filename2
 748
 749   call get_grid_size_for_one_tile(component, tile, nlon, nlat)
 750   if (size(glon(:))/=nlon) &
 751        call mpp_error (FATAL,  module_name// &
 752                        & '/get_grid_cell_centers_1D: Size of argument "glon" is not consistent with the grid size')
 753   if (size(glat(:))/=nlat) &
 754        call mpp_error (FATAL,  module_name// &
 755                        & '/get_grid_cell_centers_1D: Size of argument "glat" is not consistent with the grid size')
 756   if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
 757      call mpp_error(FATAL, module_name//'/get_grid_cell_centers_1D: Illegal component name "'// &
 758                     & trim(component)//'": must be one of ATM, LND, or OCN')
 759   endif
 760
 761   select case(get_grid_version())
 762   case(VERSION_0)
 763      select case(trim(component))
 764      case('ATM','LND')
 765         call read_data(grid_file, 'xt'//lowercase(component(1:1)), glon, no_domain=.true.)
 766         call read_data(grid_file, 'yt'//lowercase(component(1:1)), glat, no_domain=.true.)
 767      case('OCN')
 768         call read_data(grid_file, "gridlon_t", glon, no_domain=.true.)
 769         call read_data(grid_file, "gridlat_t", glat, no_domain=.true.)
 770      end select
 771   case(VERSION_1)
 772      select case(trim(component))
 773      case('ATM','LND')
 774         call read_data(grid_file, 'xt'//lowercase(component(1:1)), glon, no_domain=.true.)
 775         call read_data(grid_file, 'yt'//lowercase(component(1:1)), glat, no_domain=.true.)
 776      case('OCN')
 777         call read_data(grid_file, "grid_x_T", glon, no_domain=.true.)
 778         call read_data(grid_file, "grid_y_T", glat, no_domain=.true.)
 779      end select
 780   case(VERSION_2)
 781      ! get the name of the mosaic file for the component
 782      call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
 783      filename1=grid_dir//trim(filename1)
 784      ! get the name of the grid file for the component and tile
 785      call read_data(filename1, 'gridfiles', filename2, level=tile)
 786      filename2 = grid_dir//trim(filename2)
 787
 788      start = 1; nread = 1
 789      nread(1) = 2*nlon+1; start(2) = 2
 790      allocate( tmp(2*nlon+1,1) )
 791      call read_data(filename2, "x", tmp, start, nread, no_domain=.TRUE.)
 792      glon(1:nlon) = tmp(2:2*nlon:2,1)
 793      deallocate(tmp)
 794      allocate(tmp(1, 2*nlat+1))
 795
 796      start = 1; nread = 1
 797      nread(2) = 2*nlat+1; start(1) = 2
 798      call read_data(filename2, "y", tmp, start, nread, no_domain=.TRUE.)
 799      glat(1:nlat) = tmp(1,2:2*nlat:2)
 800      deallocate(tmp)
 801   end select
 802
 803
 804 end subroutine get_grid_cell_centers_1D
 805
 806 ! ============================================================================
 807 ! ============================================================================
 808 !> Returns grid cell centers for specified model component and mosaic tile number
 809 subroutine get_grid_cell_centers_2D(component, tile, lon, lat, domain)
 810   character(len=*), intent(in) :: component
 811   integer, intent(in) :: tile
 812   real, intent(inout) :: lon(:,:),lat(:,:)
 813   type(domain2d), intent(in), optional :: domain
 814   ! local vars
 815   character(len=MAX_FILE) :: filename1, filename2
 816   integer :: nlon, nlat
 817   integer :: i,j
 818   real, allocatable :: buffer(:),tmp(:,:)
 819   integer :: is,ie,js,je ! boundaries of our domain
 820   integer :: i0,j0 ! offsets for coordinates
 821   integer :: isg, jsg
 822   integer :: start(4), nread(4)
 823
 824   call get_grid_size_for_one_tile(component, tile, nlon, nlat)
 825   if (present(domain)) then
 826     call mpp_get_compute_domain(domain,is,ie,js,je)
 827   else
 828     is = 1 ; ie = nlon
 829     js = 1 ; je = nlat
 830     !--- domain normally should be present
 831     call mpp_error (NOTE, module_name//'/get_grid_cell_centers: domain is not present, global data will be read')
 832   endif
 833   i0 = -is+1; j0 = -js+1
 834
 835   ! verify that lon and lat sizes are consistent with the size of domain
 836   if (size(lon,1)/=ie-is+1.or.size(lon,2)/=je-js+1) &
 837        call mpp_error (FATAL, module_name// &
 838                        & '/get_grid_cell_centers: Size of array "lon" is not consistent with the domain size')
 839   if (size(lat,1)/=ie-is+1.or.size(lat,2)/=je-js+1) &
 840        call mpp_error (FATAL, module_name// &
 841                        & '/get_grid_cell_centers: Size of array "lat" is not consistent with the domain size')
 842   if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
 843      call mpp_error(FATAL, module_name//'/get_grid_cell_vertices: Illegal component name "'// &
 844                     & trim(component)//'": must be one of ATM, LND, or OCN')
 845   endif
 846
 847   select case(get_grid_version())
 848   case(VERSION_0)
 849      select case (trim(component))
 850      case('ATM','LND')
 851         allocate(buffer(max(nlon,nlat)))
 852         ! read coordinates of grid cell vertices
 853         call read_data(grid_file, 'xt'//lowercase(component(1:1)), buffer(1:nlon), no_domain=.true.)
 854         do j = js,je
 855         do i = is,ie
 856            lon(i+i0,j+j0) = buffer(i)
 857         enddo
 858         enddo
 859         call read_data(grid_file, 'yt'//lowercase(component(1:1)), buffer(1:nlat), no_domain=.true.)
 860         do j = js,je
 861         do i = is,ie
 862            lat(i+i0,j+j0) = buffer(j)
 863         enddo
 864         enddo
 865         deallocate(buffer)
 866      case('OCN')
 867         call read_data(grid_file, 'geolon_t', lon, no_domain=.not.present(domain), domain=domain )
 868         call read_data(grid_file, 'geolat_t', lat, no_domain=.not.present(domain), domain=domain )
 869      end select
 870   case(VERSION_1)
 871      select case(trim(component))
 872      case('ATM','LND')
 873         allocate(buffer(max(nlon,nlat)))
 874         ! read coordinates of grid cell vertices
 875         call read_data(grid_file, 'xt'//lowercase(component(1:1)), buffer(1:nlon), no_domain=.true.)
 876         do j = js,je
 877         do i = is,ie
 878            lon(i+i0,j+j0) = buffer(i)
 879         enddo
 880         enddo
 881         call read_data(grid_file, 'yt'//lowercase(component(1:1)), buffer(1:nlat), no_domain=.true.)
 882         do j = js,je
 883         do i = is,ie
 884            lat(i+i0,j+j0) = buffer(j)
 885         enddo
 886         enddo
 887         deallocate(buffer)
 888      case('OCN')
 889         call read_data(grid_file, 'x_T', lon, no_domain=.not.present(domain), domain=domain )
 890         call read_data(grid_file, 'y_T', lat, no_domain=.not.present(domain), domain=domain )
 891      end select
 892   case(VERSION_2) ! mosaic grid file
 893      ! get the name of the mosaic file for the component
 894      call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
 895      filename1=grid_dir//trim(filename1)
 896      ! get the name of the grid file for the component and tile
 897      call read_data(filename1, 'gridfiles', filename2, level=tile)
 898      filename2 = grid_dir//trim(filename2)
 899      if(PRESENT(domain)) then
 900         call mpp_get_global_domain(domain, xbegin=isg, ybegin=jsg)
 901         start = 1; nread = 1
 902         start(1) = 2*(is-isg+1) - 1; nread(1) = 2*(ie-is)+3
 903         start(2) = 2*(js-jsg+1) - 1; nread(2) = 2*(je-js)+3
 904         allocate(tmp(nread(1), nread(2)))
 905         call read_data(filename2, 'x', tmp, start, nread, no_domain=.TRUE.)
 906         do j = 1, je-js+1
 907            do i = 1, ie-is+1
 908               lon(i,j) = tmp(2*i,2*j)
 909            enddo
 910         enddo
 911         call read_data(filename2, 'y', tmp, start, nread, no_domain=.TRUE.)
 912         do j = 1, je-js+1
 913            do i = 1, ie-is+1
 914               lat(i,j) = tmp(2*i,2*j)
 915            enddo
 916         enddo
 917      else
 918         allocate(tmp(2*nlon+1,2*nlat+1))
 919         call read_data(filename2, 'x', tmp, no_domain=.TRUE.)
 920         do j = js,je
 921            do i = is,ie
 922               lon(i+i0,j+j0) = tmp(2*i,2*j)
 923            end do
 924         end do
 925         call read_data(filename2, 'y', tmp, no_domain=.TRUE.)
 926         do j = js,je
 927            do i = is,ie
 928               lat(i+i0,j+j0) = tmp(2*i,2*j)
 929            end do
 930         end do
 931         deallocate(tmp)
 932      endif
 933   end select
 934
 935 end subroutine get_grid_cell_centers_2D
 936
 937 subroutine get_grid_cell_centers_UG(component, tile, lon, lat, SG_domain, UG_domain)
 938   character(len=*), intent(in) :: component
 939   integer, intent(in) :: tile
 940   real, intent(inout) :: lon(:),lat(:)
 941   type(domain2d)  ,   intent(in) :: SG_domain
 942   type(domainUG)  ,   intent(in) :: UG_domain
 943   integer :: is, ie, js, je
 944   real, allocatable :: SG_lon(:,:), SG_lat(:,:)
 945
 946   call mpp_get_compute_domain(SG_domain, is, ie, js, je)
 947   allocate(SG_lon(is:ie, js:je))
 948   allocate(SG_lat(is:ie, js:je))
 949   call get_grid_cell_centers_2D(component, tile, SG_lon, SG_lat, SG_domain)
 950   call mpp_pass_SG_to_UG(UG_domain, SG_lon, lon)
 951   call mpp_pass_SG_to_UG(UG_domain, SG_lat, lat)
 952   deallocate(SG_lon, SG_lat)
 953
 954 end subroutine get_grid_cell_centers_UG
 955
 956 ! ============================================================================
 957 ! ============================================================================
 958 ! this subroutine probably does not belong in the grid_mod
 959 !> Given a model component, a layout, and (optionally) a halo size, returns a
 960 !! domain for current processor
 961 subroutine define_cube_mosaic ( component, domain, layout, halo, maskmap )
 962   character(len=*) , intent(in)    :: component
 963   type(domain2d)   , intent(inout) :: domain
 964   integer          , intent(in)    :: layout(2)
 965   integer, optional, intent(in)    :: halo
 966   logical, optional, intent(in)    :: maskmap(:,:,:)
 967
 968   ! ---- local constants
 969
 970   ! ---- local vars
 971   character(len=MAX_NAME) :: varname
 972   character(len=MAX_FILE + len(grid_dir)) :: mosaic_file
 973   integer :: ntiles     ! number of tiles
 974   integer :: ncontacts  ! number of contacts between mosaic tiles
 975   integer :: n
 976   integer :: ng, pe_pos, npes         ! halo size
 977   integer, allocatable :: nlon(:), nlat(:), global_indices(:,:)
 978   integer, allocatable :: pe_start(:), pe_end(:), layout_2d(:,:)
 979   integer, allocatable :: tile1(:),tile2(:)
 980   integer, allocatable :: is1(:),ie1(:),js1(:),je1(:)
 981   integer, allocatable :: is2(:),ie2(:),js2(:),je2(:)
 982
 983   call get_grid_ntiles(component,ntiles)
 984   allocate(nlon(ntiles), nlat(ntiles))
 985   allocate(global_indices(4,ntiles))
 986   allocate(pe_start(ntiles),pe_end(ntiles))
 987   allocate(layout_2d(2,ntiles))
 988   call get_grid_size(component,nlon,nlat)
 989
 990   pe_pos = mpp_root_pe()
 991   do n = 1, ntiles
 992      global_indices(:,n) = (/ 1, nlon(n), 1, nlat(n) /)
 993      layout_2d     (:,n) = layout
 994      if(present(maskmap)) then
 995         npes = count(maskmap(:,:,n))
 996      else
 997         npes = layout(1)*layout(2)
 998      endif
 999      pe_start(n) = pe_pos
1000      pe_end  (n) = pe_pos + npes - 1
1001      pe_pos      = pe_end(n) + 1
1002   enddo
1003
1004   varname=trim(lowercase(component))//'_mosaic_file'
1005   call read_data(grid_file,varname,mosaic_file(1:MAX_FILE))
1006   mosaic_file = grid_dir//mosaic_file(1:MAX_FILE)
1007
1008   ! get the contact information from mosaic file
1009   ncontacts = get_mosaic_ncontacts(mosaic_file)
1010   allocate(tile1(ncontacts),tile2(ncontacts))
1011   allocate(is1(ncontacts),ie1(ncontacts),js1(ncontacts),je1(ncontacts))
1012   allocate(is2(ncontacts),ie2(ncontacts),js2(ncontacts),je2(ncontacts))
1013   call get_mosaic_contact(mosaic_file, tile1, tile2, &
1014        is1, ie1, js1, je1, is2, ie2, js2, je2)
1015
1016   ng = 0
1017   if(present(halo)) ng = halo
1018   ! create the domain2d variable
1019   call mpp_define_mosaic ( global_indices, layout_2d, domain, &
1020        ntiles, ncontacts, tile1, tile2,                  &
1021        is1, ie1, js1, je1, &
1022        is2, ie2, js2, je2, &
1023        pe_start=pe_start, pe_end=pe_end, symmetry=.true.,  &
1024        shalo = ng, nhalo = ng, whalo = ng, ehalo = ng,     &
1025        maskmap = maskmap,                                  &
1026        name = trim(component)//'Cubic-Sphere Grid' )
1027
1028   deallocate(nlon,nlat,global_indices,pe_start,pe_end,layout_2d)
1029   deallocate(tile1,tile2)
1030   deallocate(is1,ie1,js1,je1)
1031   deallocate(is2,ie2,js2,je2)
1032
1033 end subroutine define_cube_mosaic
1034 #endif
1035 end module grid_mod
1036 !> @}
1037 ! close documentation grouping