!***********************************************************************
!* 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/>.
!***********************************************************************
module grid_mod
use mpp_mod, only : mpp_root_pe
use constants_mod, only : PI, radius
use fms_mod, only : uppercase, lowercase, field_exist, field_size, read_data, &
error_mesg, string, FATAL, NOTE
use fms_io_mod, only : get_great_circle_algorithm, get_global_att_value
use mosaic_mod, only : get_mosaic_ntiles, get_mosaic_xgrid_size, get_mosaic_grid_sizes, &
get_mosaic_xgrid, calc_mosaic_grid_area, calc_mosaic_grid_great_circle_area
! the following two use statement are only needed for define_cube_mosaic
use mpp_domains_mod, only : domain2d, mpp_define_mosaic, mpp_get_compute_domain, &
mpp_get_global_domain, domainUG, mpp_pass_SG_to_UG
use mosaic_mod, only : get_mosaic_ncontacts, get_mosaic_contact
implicit none;private
! ==== public interfaces =====================================================
! grid dimension inquiry subroutines
public :: get_grid_ntiles ! returns number of tiles
public :: get_grid_size ! returns horizontal sizes of the grid
! grid geometry inquiry subroutines
public :: get_grid_cell_centers
public :: get_grid_cell_vertices
! grid area inquiry subroutines
public :: get_grid_cell_area
public :: get_grid_comp_area
! decompose cubed sphere domains -- probably does not belong here, but it should
! be in some place available for component models
public :: define_cube_mosaic
! ==== end of public interfaces ==============================================
interface get_grid_size
module procedure get_grid_size_for_all_tiles
module procedure get_grid_size_for_one_tile
end interface
interface get_grid_cell_vertices
module procedure get_grid_cell_vertices_1D
module procedure get_grid_cell_vertices_2D
module procedure get_grid_cell_vertices_UG
end interface
interface get_grid_cell_centers
module procedure get_grid_cell_centers_1D
module procedure get_grid_cell_centers_2D
module procedure get_grid_cell_centers_UG
end interface
interface get_grid_cell_area
module procedure get_grid_cell_area_SG
module procedure get_grid_cell_area_UG
end interface get_grid_cell_area
interface get_grid_comp_area
module procedure get_grid_comp_area_SG
module procedure get_grid_comp_area_UG
end interface get_grid_comp_area
! ==== module constants ======================================================
character(len=*), parameter :: &
module_name = 'grid_mod'
! Include variable "version" to be written to log file.
#include<file_version.h>
character(len=*), parameter :: &
grid_dir = 'INPUT/', & ! root directory for all grid files
grid_file = 'INPUT/grid_spec.nc' ! name of the grid spec file
integer, parameter :: &
MAX_NAME = 256, & ! max length of the variable names
MAX_FILE = 1024, & ! max length of the file names
VERSION_0 = 0, &
VERSION_1 = 1, &
VERSION_2 = 2
integer, parameter :: BUFSIZE = 1048576 ! This is used to control memory usage in get_grid_comp_area
! We may change this to a namelist variable is needed.
! ==== module variables ======================================================
integer :: grid_version = -1
logical :: great_circle_algorithm = .FALSE.
logical :: first_call = .TRUE.
contains
function get_grid_version()
integer :: get_grid_version
if(first_call) then
great_circle_algorithm = get_great_circle_algorithm()
first_call = .FALSE.
endif
if(grid_version<0) then
if(field_exist(grid_file, 'geolon_t')) then
grid_version = VERSION_0
else if(field_exist(grid_file, 'x_T')) then
grid_version = VERSION_1
else if(field_exist(grid_file, 'ocn_mosaic_file') ) then
grid_version = VERSION_2
else
call error_mesg(module_name//'/get_grid_version',&
'Can''t determine the version of the grid spec: none of "x_T", "geolon_t", or "ocn_mosaic_file" exist in file "'//trim(grid_file)//'"', &
FATAL )
endif
endif
get_grid_version = grid_version
end function get_grid_version
! ============================================================================
! returns number of tiles for a given component
! ============================================================================
subroutine get_grid_ntiles(component,ntiles)
character(len=*) :: component
integer, intent(out) :: ntiles
! local vars
character(len=MAX_FILE) :: component_mosaic
select case (get_grid_version())
case(VERSION_0,VERSION_1)
ntiles = 1
case(VERSION_2)
call read_data(grid_file,trim(lowercase(component))//'_mosaic_file',component_mosaic)
ntiles = get_mosaic_ntiles(grid_dir//trim(component_mosaic))
end select
end subroutine get_grid_ntiles
! ============================================================================
! returns size of the grid for each of the tiles
! ============================================================================
subroutine get_grid_size_for_all_tiles(component,nx,ny)
character(len=*) :: component
integer, intent(inout) :: nx(:),ny(:)
! local vars
integer :: siz(4) ! for the size of external fields
character(len=MAX_NAME) :: varname1, varname2
character(len=MAX_FILE) :: component_mosaic
varname1 = 'AREA_'//trim(uppercase(component))
varname2 = trim(lowercase(component))//'_mosaic_file'
select case (get_grid_version())
case(VERSION_0,VERSION_1)
call field_size(grid_file, varname1, siz)
nx(1) = siz(1); ny(1)=siz(2)
case(VERSION_2) ! mosaic file
call read_data(grid_file,varname2, component_mosaic)
call get_mosaic_grid_sizes(grid_dir//trim(component_mosaic),nx,ny)
end select
end subroutine get_grid_size_for_all_tiles
! ============================================================================
! returns size of the grid for one of the tiles
! ============================================================================
subroutine get_grid_size_for_one_tile(component,tile,nx,ny)
character(len=*) :: component
integer, intent(in) :: tile
integer, intent(inout) :: nx,ny
! local vars
integer, allocatable :: nnx(:), nny(:)
integer :: ntiles
call get_grid_ntiles(component, ntiles)
if(tile>0.and.tile<=ntiles) then
allocate(nnx(ntiles),nny(ntiles))
call get_grid_size_for_all_tiles(component,nnx,nny)
nx = nnx(tile); ny = nny(tile)
deallocate(nnx,nny)
else
call error_mesg('get_grid_size',&
'requested tile index '//trim(string(tile))//' is out of bounds (1:'//trim(string(ntiles))//')',&
FATAL)
endif
end subroutine get_grid_size_for_one_tile
! ============================================================================
! return grid cell area for the specified model component and tile
! ============================================================================
subroutine get_grid_cell_area_SG(component, tile, cellarea, domain)
character(len=*), intent(in) :: component
integer , intent(in) :: tile
real , intent(inout) :: cellarea(:,:)
type(domain2d) , intent(in), optional :: domain
! local vars
integer :: nlon, nlat
real, allocatable :: glonb(:,:), glatb(:,:)
select case(get_grid_version())
case(VERSION_0,VERSION_1)
select case(trim(component))
case('LND')
call read_data(grid_file, 'AREA_LND_CELL', cellarea, &
no_domain=.not.present(domain), domain=domain)
case('ATM','OCN')
call read_data(grid_file, 'AREA_'//trim(uppercase(component)),cellarea,&
no_domain=.not.present(domain),domain=domain)
case default
call error_mesg(module_name//'/get_grid_cell_area',&
'Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN',&
FATAL)
end select
! convert area to m2
cellarea = cellarea*4.*PI*radius**2
case(VERSION_2)
if (present(domain)) then
call mpp_get_compute_domain(domain,xsize=nlon,ysize=nlat)
else
call get_grid_size(component,tile,nlon,nlat)
endif
allocate(glonb(nlon+1,nlat+1),glatb(nlon+1,nlat+1))
call get_grid_cell_vertices(component, tile, glonb, glatb, domain)
if (great_circle_algorithm) then
call calc_mosaic_grid_great_circle_area(glonb*pi/180.0, glatb*pi/180.0, cellarea)
else
call calc_mosaic_grid_area(glonb*pi/180.0, glatb*pi/180.0, cellarea)
end if
deallocate(glonb,glatb)
end select
end subroutine get_grid_cell_area_SG
! ============================================================================
! get the area of the component per grid cell
! ============================================================================
subroutine get_grid_comp_area_SG(component,tile,area,domain)
character(len=*) :: component
integer, intent(in) :: tile
real, intent(inout) :: area(:,:)
type(domain2d), intent(in), optional :: domain
! local vars
integer :: n_xgrid_files ! number of exchange grid files in the mosaic
integer :: siz(4), nxgrid
integer :: i,j,m,n
integer, allocatable :: i1(:), j1(:), i2(:), j2(:)
real, allocatable :: xgrid_area(:)
real, allocatable :: rmask(:,:)
character(len=MAX_NAME) :: &
xgrid_name, & ! name of the variable holding xgrid names
tile_name, & ! name of the tile
xgrid_file, & ! name of the current xgrid file
mosaic_name,& ! name of the mosaic
mosaic_file,&
tilefile
character(len=4096) :: attvalue
character(len=MAX_NAME), allocatable :: nest_tile_name(:)
character(len=MAX_NAME) :: varname1, varname2
integer :: is,ie,js,je ! boundaries of our domain
integer :: i0, j0 ! offsets for x and y, respectively
integer :: num_nest_tile, ntiles
logical :: is_nest
integer :: found_xgrid_files ! how many xgrid files we actually found in the grid spec
integer :: ibegin, iend, bsize, l
select case (get_grid_version())
case(VERSION_0,VERSION_1)
select case(component)
case('ATM')
call read_data(grid_file,'AREA_ATM',area, no_domain=.not.present(domain),domain=domain)
case('OCN')
allocate(rmask(size(area,1),size(area,2)))
call read_data(grid_file,'AREA_OCN',area, no_domain=.not.present(domain),domain=domain)
call read_data(grid_file,'wet', rmask,no_domain=.not.present(domain),domain=domain)
area = area*rmask
deallocate(rmask)
case('LND')
call read_data(grid_file,'AREA_LND',area,no_domain=.not.present(domain),domain=domain)
case default
call error_mesg(module_name//'/get_grid_comp_area',&
'Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN',&
FATAL)
end select
case(VERSION_2) ! mosaic gridspec
select case (component)
case ('ATM')
! just read the grid cell area and return
call get_grid_cell_area(component,tile,area)
return
case ('LND')
xgrid_name = 'aXl_file'
call read_data(grid_file, 'lnd_mosaic', mosaic_name)
tile_name = trim(mosaic_name)//'_tile'//char(tile+ichar('0'))
case ('OCN')
xgrid_name = 'aXo_file'
call read_data(grid_file, 'ocn_mosaic', mosaic_name)
tile_name = trim(mosaic_name)//'_tile'//char(tile+ichar('0'))
case default
call error_mesg(module_name//'/get_grid_comp_area',&
'Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN',&
FATAL)
end select
! get the boundaries of the requested domain
if(present(domain)) then
call mpp_get_compute_domain(domain,is,ie,js,je)
i0 = 1-is ; j0=1-js
else
call get_grid_size(component,tile,ie,je)
is = 1 ; i0 = 0
js = 1 ; j0 = 0
endif
if (size(area,1)/=ie-is+1.or.size(area,2)/=je-js+1) &
call error_mesg(module_name//'/get_grid_comp_area',&
'size of the output argument "area" is not consistent with the domain',FATAL)
! find the nest tile
call read_data(grid_file, 'atm_mosaic', mosaic_name)
call read_data(grid_file,'atm_mosaic_file',mosaic_file)
mosaic_file = grid_dir//trim(mosaic_file)
ntiles = get_mosaic_ntiles(trim(mosaic_file))
allocate(nest_tile_name(ntiles))
num_nest_tile = 0
do n = 1, ntiles
call read_data(mosaic_file, 'gridfiles', tilefile, level=n)
tilefile = grid_dir//trim(tilefile)
if( get_global_att_value(tilefile, "nest_grid", attvalue) ) then
if(trim(attvalue) == "TRUE") then
num_nest_tile = num_nest_tile + 1
nest_tile_name(num_nest_tile) = trim(mosaic_name)//'_tile'//char(n+ichar('0'))
else if(trim(attvalue) .NE. "FALSE") then
call error_mesg(module_name//'/get_grid_comp_area', 'value of global attribute nest_grid in file'// &
trim(tilefile)//' should be TRUE of FALSE', FATAL)
endif
end if
end do
area(:,:) = 0.
if(field_exist(grid_file,xgrid_name)) then
! get the number of the exchange-grid files
call field_size(grid_file,xgrid_name,siz)
n_xgrid_files = siz(2)
found_xgrid_files = 0
! loop through all exchange grid files
do n = 1, n_xgrid_files
! get the name of the current exchange grid file
call read_data(grid_file,xgrid_name,xgrid_file,level=n)
! skip the rest of the loop if the name of the current tile isn't found
! in the file name, but check this only if there is more than 1 tile
if(n_xgrid_files>1) then
if(index(xgrid_file,trim(tile_name))==0) cycle
endif
found_xgrid_files = found_xgrid_files + 1
!---make sure the atmosphere grid is not a nested grid
is_nest = .false.
do m = 1, num_nest_tile
if(index(xgrid_file, trim(nest_tile_name(m))) .NE. 0) then
is_nest = .true.
exit
end if
end do
if(is_nest) cycle
! finally read the exchange grid
nxgrid = get_mosaic_xgrid_size(grid_dir//xgrid_file)
if(nxgrid < BUFSIZE) then
allocate(i1(nxgrid), j1(nxgrid), i2(nxgrid), j2(nxgrid), xgrid_area(nxgrid))
else
allocate(i1(BUFSIZE), j1(BUFSIZE), i2(BUFSIZE), j2(BUFSIZE), xgrid_area(BUFSIZE))
endif
ibegin = 1
do l = 1,nxgrid,BUFSIZE
bsize = min(BUFSIZE, nxgrid-l+1)
iend = ibegin + bsize - 1
call get_mosaic_xgrid(grid_dir//xgrid_file, i1(1:bsize), j1(1:bsize), i2(1:bsize), j2(1:bsize), &
xgrid_area(1:bsize), ibegin, iend)
! and sum the exchange grid areas
do m = 1, bsize
i = i2(m); j = j2(m)
if (i<is.or.i>ie) cycle
if (j<js.or.j>je) cycle
area(i+i0,j+j0) = area(i+i0,j+j0) + xgrid_area(m)
end do
ibegin = iend + 1
enddo
deallocate(i1, j1, i2, j2, xgrid_area)
enddo
if (found_xgrid_files == 0) &
call error_mesg('get_grid_comp_area', 'no xgrid files were found for component '&
//trim(component)//' (mosaic name is '//trim(mosaic_name)//')', FATAL)
endif
deallocate(nest_tile_name)
end select ! version
! convert area to m2
area = area*4.*PI*radius**2
end subroutine get_grid_comp_area_SG
!======================================================================
subroutine get_grid_cell_area_UG(component, tile, cellarea, SG_domain, UG_domain)
character(len=*), intent(in) :: component
integer , intent(in) :: tile
real , intent(inout) :: cellarea(:)
type(domain2d) , intent(in) :: SG_domain
type(domainUG) , intent(in) :: UG_domain
integer :: is, ie, js, je
real, allocatable :: SG_area(:,:)
call mpp_get_compute_domain(SG_domain, is, ie, js, je)
allocate(SG_area(is:ie, js:je))
call get_grid_cell_area_SG(component, tile, SG_area, SG_domain)
call mpp_pass_SG_to_UG(UG_domain, SG_area, cellarea)
deallocate(SG_area)
end subroutine get_grid_cell_area_UG
subroutine get_grid_comp_area_UG(component, tile, area, SG_domain, UG_domain)
character(len=*), intent(in) :: component
integer , intent(in) :: tile
real , intent(inout) :: area(:)
type(domain2d) , intent(in) :: SG_domain
type(domainUG) , intent(in) :: UG_domain
integer :: is, ie, js, je
real, allocatable :: SG_area(:,:)
call mpp_get_compute_domain(SG_domain, is, ie, js, je)
allocate(SG_area(is:ie, js:je))
call get_grid_comp_area_SG(component, tile, SG_area, SG_domain)
call mpp_pass_SG_to_UG(UG_domain, SG_area, area)
deallocate(SG_area)
end subroutine get_grid_comp_area_UG
! ============================================================================
! returns arrays of global grid cell boundaries for given model component and
! mosaic tile number.
! NOTE that in case of non-lat-lon grid the returned coordinates may have be not so
! meaningful, by the very nature of such grids. But presumably these 1D coordinate
! arrays are good enough for diag axis and such.
! ============================================================================
subroutine get_grid_cell_vertices_1D(component, tile, glonb, glatb)
character(len=*), intent(in) :: component
integer, intent(in) :: tile
real, intent(inout) :: glonb(:),glatb(:)
integer :: nlon, nlat
integer :: start(4), nread(4)
real, allocatable :: tmp(:,:), x_vert_t(:,:,:), y_vert_t(:,:,:)
character(len=MAX_FILE) :: filename1, filename2
call get_grid_size_for_one_tile(component, tile, nlon, nlat)
if (size(glonb(:))/=nlon+1) &
call error_mesg ( module_name//'/get_grid_cell_vertices_1D',&
'Size of argument "glonb" is not consistent with the grid size',FATAL)
if (size(glatb(:))/=nlat+1) &
call error_mesg ( module_name//'/get_grid_cell_vertices_1D',&
'Size of argument "glatb" is not consistent with the grid size',FATAL)
if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
call error_mesg(module_name//'/get_grid_cell_vertices_1D',&
'Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN',&
FATAL)
endif
select case(get_grid_version())
case(VERSION_0)
select case(trim(component))
case('ATM','LND')
call read_data(grid_file, 'xb'//lowercase(component(1:1)), glonb, no_domain=.true.)
call read_data(grid_file, 'yb'//lowercase(component(1:1)), glatb, no_domain=.true.)
case('OCN')
call read_data(grid_file, "gridlon_vert_t", glonb, no_domain=.true.)
call read_data(grid_file, "gridlat_vert_t", glatb, no_domain=.true.)
end select
case(VERSION_1)
select case(trim(component))
case('ATM','LND')
call read_data(grid_file, 'xb'//lowercase(component(1:1)), glonb, no_domain=.true.)
call read_data(grid_file, 'yb'//lowercase(component(1:1)), glatb, no_domain=.true.)
case('OCN')
allocate (x_vert_t(nlon,1,2), y_vert_t(1,nlat,2) )
start = 1; nread = 1
nread(1) = nlon; nread(2) = 1; start(3) = 1
call read_data(grid_file, "x_vert_T", x_vert_t(:,:,1), start, nread, no_domain=.TRUE.)
nread(1) = nlon; nread(2) = 1; start(3) = 2
call read_data(grid_file, "x_vert_T", x_vert_t(:,:,2), start, nread, no_domain=.TRUE.)
nread(1) = 1; nread(2) = nlat; start(3) = 1
call read_data(grid_file, "y_vert_T", y_vert_t(:,:,1), start, nread, no_domain=.TRUE.)
nread(1) = 1; nread(2) = nlat; start(3) = 4
call read_data(grid_file, "y_vert_T", y_vert_t(:,:,2), start, nread, no_domain=.TRUE.)
glonb(1:nlon) = x_vert_t(1:nlon,1,1)
glonb(nlon+1) = x_vert_t(nlon,1,2)
glatb(1:nlat) = y_vert_t(1,1:nlat,1)
glatb(nlat+1) = y_vert_t(1,nlat,2)
deallocate(x_vert_t, y_vert_t)
end select
case(VERSION_2)
! get the name of the mosaic file for the component
call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
filename1=grid_dir//trim(filename1)
! get the name of the grid file for the component and tile
call read_data(filename1, 'gridfiles', filename2, level=tile)
filename2 = grid_dir//trim(filename2)
start = 1; nread = 1
nread(1) = 2*nlon+1
allocate( tmp(2*nlon+1,1) )
call read_data(filename2, "x", tmp, start, nread, no_domain=.TRUE.)
glonb(1:nlon+1) = tmp(1:2*nlon+1:2,1)
deallocate(tmp)
allocate(tmp(1,2*nlat+1))
start = 1; nread = 1
nread(2) = 2*nlat+1
call read_data(filename2, "y", tmp, start, nread, no_domain=.TRUE.)
glatb(1:nlat+1) = tmp(1,1:2*nlat+1:2)
deallocate(tmp)
end select
end subroutine get_grid_cell_vertices_1D
! ============================================================================
! returns cell vertices for the specified model component and mosaic tile number
! ============================================================================
subroutine get_grid_cell_vertices_2D(component, tile, lonb, latb, domain)
character(len=*), intent(in) :: component
integer, intent(in) :: tile
real, intent(inout) :: lonb(:,:),latb(:,:)
type(domain2d), optional, intent(in) :: domain
! local vars
character(len=MAX_FILE) :: filename1, filename2
integer :: nlon, nlat
integer :: i,j
real, allocatable :: buffer(:), tmp(:,:), x_vert_t(:,:,:), y_vert_t(:,:,:)
integer :: is,ie,js,je ! boundaries of our domain
integer :: i0,j0 ! offsets for coordinates
integer :: isg, jsg
integer :: start(4), nread(4)
call get_grid_size_for_one_tile(component, tile, nlon, nlat)
if (present(domain)) then
call mpp_get_compute_domain(domain,is,ie,js,je)
else
is = 1 ; ie = nlon
js = 1 ; je = nlat
!--- domain normally should be present
call error_mesg ( module_name//'/get_grid_cell_vertices',&
'domain is not present, global data will be read', NOTE)
endif
i0 = -is+1; j0 = -js+1
! verify that lonb and latb sizes are consistent with the size of domain
if (size(lonb,1)/=ie-is+2.or.size(lonb,2)/=je-js+2) &
call error_mesg ( module_name//'/get_grid_cell_vertices',&
'Size of argument "lonb" is not consistent with the domain size',FATAL)
if (size(latb,1)/=ie-is+2.or.size(latb,2)/=je-js+2) &
call error_mesg ( module_name//'/get_grid_cell_vertices',&
'Size of argument "latb" is not consistent with the domain size',FATAL)
if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
call error_mesg(module_name//'/get_grid_cell_vertices',&
'Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN',&
FATAL)
endif
select case(get_grid_version())
case(VERSION_0)
select case(component)
case('ATM','LND')
allocate(buffer(max(nlon,nlat)+1))
! read coordinates of grid cell vertices
call read_data(grid_file, 'xb'//lowercase(component(1:1)), buffer(1:nlon+1), no_domain=.true.)
do j = js, je+1
do i = is, ie+1
lonb(i+i0,j+j0) = buffer(i)
enddo
enddo
call read_data(grid_file, 'yb'//lowercase(component(1:1)), buffer(1:nlat+1), no_domain=.true.)
do j = js, je+1
do i = is, ie+1
latb(i+i0,j+j0) = buffer(j)
enddo
enddo
deallocate(buffer)
case('OCN')
if (present(domain)) then
start = 1; nread = 1
start(1) = is; start(2) = js
nread(1) = ie-is+2; nread(2) = je-js+2
call read_data(grid_file, 'geolon_vert_t', lonb, start, nread, no_domain=.true. )
call read_data(grid_file, 'geolat_vert_t', latb, start, nread, no_domain=.true. )
else
call read_data(grid_file, 'geolon_vert_t', lonb, no_domain=.TRUE. )
call read_data(grid_file, 'geolat_vert_t', latb, no_domain=.TRUE. )
endif
end select
case(VERSION_1)
select case(component)
case('ATM','LND')
allocate(buffer(max(nlon,nlat)+1))
! read coordinates of grid cell vertices
call read_data(grid_file, 'xb'//lowercase(component(1:1)), buffer(1:nlon+1), no_domain=.true.)
do j = js, je+1
do i = is, ie+1
lonb(i+i0,j+j0) = buffer(i)
enddo
enddo
call read_data(grid_file, 'yb'//lowercase(component(1:1)), buffer(1:nlat+1), no_domain=.true.)
do j = js, je+1
do i = is, ie+1
latb(i+i0,j+j0) = buffer(j)
enddo
enddo
deallocate(buffer)
case('OCN')
nlon=ie-is+1; nlat=je-js+1
allocate (x_vert_t(nlon,nlat,4), y_vert_t(nlon,nlat,4) )
call read_data(grid_file, 'x_vert_T', x_vert_t, no_domain=.not.present(domain), domain=domain )
call read_data(grid_file, 'y_vert_T', y_vert_t, no_domain=.not.present(domain), domain=domain )
lonb(1:nlon,1:nlat) = x_vert_t(1:nlon,1:nlat,1)
lonb(nlon+1,1:nlat) = x_vert_t(nlon,1:nlat,2)
lonb(1:nlon,nlat+1) = x_vert_t(1:nlon,nlat,4)
lonb(nlon+1,nlat+1) = x_vert_t(nlon,nlat,3)
latb(1:nlon,1:nlat) = y_vert_t(1:nlon,1:nlat,1)
latb(nlon+1,1:nlat) = y_vert_t(nlon,1:nlat,2)
latb(1:nlon,nlat+1) = y_vert_t(1:nlon,nlat,4)
latb(nlon+1,nlat+1) = y_vert_t(nlon,nlat,3)
deallocate(x_vert_t, y_vert_t)
end select
case(VERSION_2)
! get the name of the mosaic file for the component
call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
filename1=grid_dir//trim(filename1)
! get the name of the grid file for the component and tile
call read_data(filename1, 'gridfiles', filename2, level=tile)
filename2 = grid_dir//trim(filename2)
if(PRESENT(domain)) then
call mpp_get_global_domain(domain, xbegin=isg, ybegin=jsg)
start = 1; nread = 1
start(1) = 2*(is-isg+1) - 1; nread(1) = 2*(ie-is)+3
start(2) = 2*(js-jsg+1) - 1; nread(2) = 2*(je-js)+3
allocate(tmp(nread(1), nread(2)) )
call read_data(filename2, 'x', tmp, start, nread, no_domain=.TRUE.)
do j = 1, je-js+2
do i = 1, ie-is+2
lonb(i,j) = tmp(2*i-1,2*j-1)
enddo
enddo
call read_data(filename2, 'y', tmp, start, nread, no_domain=.TRUE.)
do j = 1, je-js+2
do i = 1, ie-is+2
latb(i,j) = tmp(2*i-1,2*j-1)
enddo
enddo
else
allocate(tmp(2*nlon+1,2*nlat+1))
call read_data(filename2, 'x', tmp, no_domain=.TRUE.)
do j = js, je+1
do i = is, ie+1
lonb(i+i0,j+j0) = tmp(2*i-1,2*j-1)
end do
end do
call read_data(filename2, 'y', tmp, no_domain=.TRUE.)
do j = js, je+1
do i = is, ie+1
latb(i+i0,j+j0) = tmp(2*i-1,2*j-1)
end do
end do
endif
deallocate(tmp)
end select
end subroutine get_grid_cell_vertices_2D
subroutine get_grid_cell_vertices_UG(component, tile, lonb, latb, SG_domain, UG_domain)
character(len=*), intent(in) :: component
integer, intent(in) :: tile
real, intent(inout) :: lonb(:,:),latb(:,:) ! The second dimension is 4
type(domain2d) , intent(in) :: SG_domain
type(domainUG) , intent(in) :: UG_domain
integer :: is, ie, js, je, i, j
real, allocatable :: SG_lonb(:,:), SG_latb(:,:), tmp(:,:,:)
call mpp_get_compute_domain(SG_domain, is, ie, js, je)
allocate(SG_lonb(is:ie+1, js:je+1))
allocate(SG_latb(is:ie+1, js:je+1))
allocate(tmp(is:ie,js:je,4))
call get_grid_cell_vertices_2D(component, tile, SG_lonb, SG_latb, SG_domain)
do j = js, je
do i = is, ie
tmp(i,j,1) = SG_lonb(i,j)
tmp(i,j,2) = SG_lonb(i+1,j)
tmp(i,j,3) = SG_lonb(i+1,j+1)
tmp(i,j,4) = SG_lonb(i,j+1)
enddo
enddo
call mpp_pass_SG_to_UG(UG_domain, tmp, lonb)
do j = js, je
do i = is, ie
tmp(i,j,1) = SG_latb(i,j)
tmp(i,j,2) = SG_latb(i+1,j)
tmp(i,j,3) = SG_latb(i+1,j+1)
tmp(i,j,4) = SG_latb(i,j+1)
enddo
enddo
call mpp_pass_SG_to_UG(UG_domain, tmp, latb)
deallocate(SG_lonb, SG_latb, tmp)
end subroutine get_grid_cell_vertices_UG
! ============================================================================
! returns global coordinate arrays fro given model component and mosaic tile number
! NOTE that in case of non-lat-lon grid those coordinates may have be not so
! meaningful, by the very nature of such grids. But presumably these 1D coordinate
! arrays are good enough for diag axis and such.
! ============================================================================
subroutine get_grid_cell_centers_1D(component, tile, glon, glat)
character(len=*), intent(in) :: component
integer, intent(in) :: tile
real, intent(inout) :: glon(:),glat(:)
integer :: nlon, nlat
integer :: start(4), nread(4)
real, allocatable :: tmp(:,:)
character(len=MAX_FILE) :: filename1, filename2
call get_grid_size_for_one_tile(component, tile, nlon, nlat)
if (size(glon(:))/=nlon) &
call error_mesg ( module_name//'/get_grid_cell_centers_1D',&
'Size of argument "glon" is not consistent with the grid size',FATAL)
if (size(glat(:))/=nlat) &
call error_mesg ( module_name//'/get_grid_cell_centers_1D',&
'Size of argument "glat" is not consistent with the grid size',FATAL)
if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
call error_mesg(module_name//'/get_grid_cell_centers_1D',&
'Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN',&
FATAL)
endif
select case(get_grid_version())
case(VERSION_0)
select case(trim(component))
case('ATM','LND')
call read_data(grid_file, 'xt'//lowercase(component(1:1)), glon, no_domain=.true.)
call read_data(grid_file, 'yt'//lowercase(component(1:1)), glat, no_domain=.true.)
case('OCN')
call read_data(grid_file, "gridlon_t", glon, no_domain=.true.)
call read_data(grid_file, "gridlat_t", glat, no_domain=.true.)
end select
case(VERSION_1)
select case(trim(component))
case('ATM','LND')
call read_data(grid_file, 'xt'//lowercase(component(1:1)), glon, no_domain=.true.)
call read_data(grid_file, 'yt'//lowercase(component(1:1)), glat, no_domain=.true.)
case('OCN')
call read_data(grid_file, "grid_x_T", glon, no_domain=.true.)
call read_data(grid_file, "grid_y_T", glat, no_domain=.true.)
end select
case(VERSION_2)
! get the name of the mosaic file for the component
call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
filename1=grid_dir//trim(filename1)
! get the name of the grid file for the component and tile
call read_data(filename1, 'gridfiles', filename2, level=tile)
filename2 = grid_dir//trim(filename2)
start = 1; nread = 1
nread(1) = 2*nlon+1; start(2) = 2
allocate( tmp(2*nlon+1,1) )
call read_data(filename2, "x", tmp, start, nread, no_domain=.TRUE.)
glon(1:nlon) = tmp(2:2*nlon:2,1)
deallocate(tmp)
allocate(tmp(1, 2*nlat+1))
start = 1; nread = 1
nread(2) = 2*nlat+1; start(1) = 2
call read_data(filename2, "y", tmp, start, nread, no_domain=.TRUE.)
glat(1:nlat) = tmp(1,2:2*nlat:2)
deallocate(tmp)
end select
end subroutine get_grid_cell_centers_1D
! ============================================================================
! returns grid cell centers for specified model component and mosaic tile number
! ============================================================================
subroutine get_grid_cell_centers_2D(component, tile, lon, lat, domain)
character(len=*), intent(in) :: component
integer, intent(in) :: tile
real, intent(inout) :: lon(:,:),lat(:,:)
type(domain2d), intent(in), optional :: domain
! local vars
character(len=MAX_NAME) :: varname
character(len=MAX_FILE) :: filename1, filename2
integer :: nlon, nlat
integer :: i,j
real, allocatable :: buffer(:),tmp(:,:)
integer :: is,ie,js,je ! boundaries of our domain
integer :: i0,j0 ! offsets for coordinates
integer :: isg, jsg
integer :: start(4), nread(4)
call get_grid_size_for_one_tile(component, tile, nlon, nlat)
if (present(domain)) then
call mpp_get_compute_domain(domain,is,ie,js,je)
else
is = 1 ; ie = nlon
js = 1 ; je = nlat
!--- domain normally should be present
call error_mesg ( module_name//'/get_grid_cell_centers',&
'domain is not present, global data will be read', NOTE)
endif
i0 = -is+1; j0 = -js+1
! verify that lon and lat sizes are consistent with the size of domain
if (size(lon,1)/=ie-is+1.or.size(lon,2)/=je-js+1) &
call error_mesg ( module_name//'/get_grid_cell_centers',&
'Size of array "lon" is not consistent with the domain size',&
FATAL )
if (size(lat,1)/=ie-is+1.or.size(lat,2)/=je-js+1) &
call error_mesg ( module_name//'/get_grid_cell_centers',&
'Size of array "lat" is not consistent with the domain size',&
FATAL )
if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
call error_mesg(module_name//'/get_grid_cell_vertices',&
'Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN',&
FATAL)
endif
select case(get_grid_version())
case(VERSION_0)
select case (trim(component))
case('ATM','LND')
allocate(buffer(max(nlon,nlat)))
! read coordinates of grid cell vertices
call read_data(grid_file, 'xt'//lowercase(component(1:1)), buffer(1:nlon), no_domain=.true.)
do j = js,je
do i = is,ie
lon(i+i0,j+j0) = buffer(i)
enddo
enddo
call read_data(grid_file, 'yt'//lowercase(component(1:1)), buffer(1:nlat), no_domain=.true.)
do j = js,je
do i = is,ie
lat(i+i0,j+j0) = buffer(j)
enddo
enddo
deallocate(buffer)
case('OCN')
call read_data(grid_file, 'geolon_t', lon, no_domain=.not.present(domain), domain=domain )
call read_data(grid_file, 'geolat_t', lat, no_domain=.not.present(domain), domain=domain )
end select
case(VERSION_1)
select case(trim(component))
case('ATM','LND')
allocate(buffer(max(nlon,nlat)))
! read coordinates of grid cell vertices
call read_data(grid_file, 'xt'//lowercase(component(1:1)), buffer(1:nlon), no_domain=.true.)
do j = js,je
do i = is,ie
lon(i+i0,j+j0) = buffer(i)
enddo
enddo
call read_data(grid_file, 'yt'//lowercase(component(1:1)), buffer(1:nlat), no_domain=.true.)
do j = js,je
do i = is,ie
lat(i+i0,j+j0) = buffer(j)
enddo
enddo
deallocate(buffer)
case('OCN')
call read_data(grid_file, 'x_T', lon, no_domain=.not.present(domain), domain=domain )
call read_data(grid_file, 'y_T', lat, no_domain=.not.present(domain), domain=domain )
end select
case(VERSION_2) ! mosaic grid file
! get the name of the mosaic file for the component
call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
filename1=grid_dir//trim(filename1)
! get the name of the grid file for the component and tile
call read_data(filename1, 'gridfiles', filename2, level=tile)
filename2 = grid_dir//trim(filename2)
if(PRESENT(domain)) then
call mpp_get_global_domain(domain, xbegin=isg, ybegin=jsg)
start = 1; nread = 1
start(1) = 2*(is-isg+1) - 1; nread(1) = 2*(ie-is)+3
start(2) = 2*(js-jsg+1) - 1; nread(2) = 2*(je-js)+3
allocate(tmp(nread(1), nread(2)))
call read_data(filename2, 'x', tmp, start, nread, no_domain=.TRUE.)
do j = 1, je-js+1
do i = 1, ie-is+1
lon(i,j) = tmp(2*i,2*j)
enddo
enddo
call read_data(filename2, 'y', tmp, start, nread, no_domain=.TRUE.)
do j = 1, je-js+1
do i = 1, ie-is+1
lat(i,j) = tmp(2*i,2*j)
enddo
enddo
else
allocate(tmp(2*nlon+1,2*nlat+1))
call read_data(filename2, 'x', tmp, no_domain=.TRUE.)
do j = js,je
do i = is,ie
lon(i+i0,j+j0) = tmp(2*i,2*j)
end do
end do
call read_data(filename2, 'y', tmp, no_domain=.TRUE.)
do j = js,je
do i = is,ie
lat(i+i0,j+j0) = tmp(2*i,2*j)
end do
end do
deallocate(tmp)
endif
end select
end subroutine get_grid_cell_centers_2D
subroutine get_grid_cell_centers_UG(component, tile, lon, lat, SG_domain, UG_domain)
character(len=*), intent(in) :: component
integer, intent(in) :: tile
real, intent(inout) :: lon(:),lat(:)
type(domain2d) , intent(in) :: SG_domain
type(domainUG) , intent(in) :: UG_domain
integer :: is, ie, js, je
real, allocatable :: SG_lon(:,:), SG_lat(:,:)
call mpp_get_compute_domain(SG_domain, is, ie, js, je)
allocate(SG_lon(is:ie, js:je))
allocate(SG_lat(is:ie, js:je))
call get_grid_cell_centers_2D(component, tile, SG_lon, SG_lat, SG_domain)
call mpp_pass_SG_to_UG(UG_domain, SG_lon, lon)
call mpp_pass_SG_to_UG(UG_domain, SG_lat, lat)
deallocate(SG_lon, SG_lat)
end subroutine get_grid_cell_centers_UG
! ============================================================================
! given a model component, a layout, and (optionally) a halo size, returns a
! domain for current processor
! ============================================================================
! this subroutine probably does not belong in the grid_mod
subroutine define_cube_mosaic ( component, domain, layout, halo, maskmap )
character(len=*) , intent(in) :: component
type(domain2d) , intent(inout) :: domain
integer , intent(in) :: layout(2)
integer, optional, intent(in) :: halo
logical, optional, intent(in) :: maskmap(:,:,:)
! ---- local constants
! ---- local vars
character(len=MAX_NAME) :: varname
character(len=MAX_FILE) :: mosaic_file
integer :: ntiles ! number of tiles
integer :: ncontacts ! number of contacts between mosaic tiles
integer :: n
integer :: ng, pe_pos, npes ! halo size
integer, allocatable :: nlon(:), nlat(:), global_indices(:,:)
integer, allocatable :: pe_start(:), pe_end(:), layout_2d(:,:)
integer, allocatable :: tile1(:),tile2(:)
integer, allocatable :: is1(:),ie1(:),js1(:),je1(:)
integer, allocatable :: is2(:),ie2(:),js2(:),je2(:)
call get_grid_ntiles(component,ntiles)
allocate(nlon(ntiles), nlat(ntiles))
allocate(global_indices(4,ntiles))
allocate(pe_start(ntiles),pe_end(ntiles))
allocate(layout_2d(2,ntiles))
call get_grid_size(component,nlon,nlat)
pe_pos = mpp_root_pe()
do n = 1, ntiles
global_indices(:,n) = (/ 1, nlon(n), 1, nlat(n) /)
layout_2d (:,n) = layout
if(present(maskmap)) then
npes = count(maskmap(:,:,n))
else
npes = layout(1)*layout(2)
endif
pe_start(n) = pe_pos
pe_end (n) = pe_pos + npes - 1
pe_pos = pe_end(n) + 1
enddo
varname=trim(lowercase(component))//'_mosaic_file'
call read_data(grid_file,varname,mosaic_file)
mosaic_file = grid_dir//mosaic_file
! get the contact information from mosaic file
ncontacts = get_mosaic_ncontacts(mosaic_file)
allocate(tile1(ncontacts),tile2(ncontacts))
allocate(is1(ncontacts),ie1(ncontacts),js1(ncontacts),je1(ncontacts))
allocate(is2(ncontacts),ie2(ncontacts),js2(ncontacts),je2(ncontacts))
call get_mosaic_contact(mosaic_file, tile1, tile2, &
is1, ie1, js1, je1, is2, ie2, js2, je2)
ng = 0
if(present(halo)) ng = halo
! create the domain2d variable
call mpp_define_mosaic ( global_indices, layout_2d, domain, &
ntiles, ncontacts, tile1, tile2, &
is1, ie1, js1, je1, &
is2, ie2, js2, je2, &
pe_start=pe_start, pe_end=pe_end, symmetry=.true., &
shalo = ng, nhalo = ng, whalo = ng, ehalo = ng, &
maskmap = maskmap, &
name = trim(component)//'Cubic-Sphere Grid' )
deallocate(nlon,nlat,global_indices,pe_start,pe_end,layout_2d)
deallocate(tile1,tile2)
deallocate(is1,ie1,js1,je1)
deallocate(is2,ie2,js2,je2)
end subroutine define_cube_mosaic
end module grid_mod