subroutine MAPL_ReadTilingNC4(File, GridName, im, jm, nx, ny, n_Grids, n_tiles, iTable, rTable, N_PfafCat, AVR,rc)
character(*), intent(IN) :: File
character(*), optional, intent(out) :: GridName(:)
integer, optional, intent(out) :: IM(:), JM(:)
integer, optional, intent(out) :: nx, ny, n_Grids, n_tiles
integer, optional, allocatable, intent(out) :: iTable(:,:)
real(kind=REAL64), optional, allocatable, intent(out) :: rTable(:,:)
integer, optional, intent(out) :: N_PfafCat
real, optional, pointer, intent(out) :: AVR(:,:) ! used by GEOSgcm
integer, optional, intent(out) :: rc
type (Attribute), pointer :: ref
character(len=:), allocatable :: attr
type (NetCDF4_FileFormatter) :: formatter
type (FileMetadata) :: meta
character(len=4) :: ocn_str
integer :: ng, ntile, status, ll
class(*), pointer :: attr_val(:)
class(*), pointer :: char_val
integer, allocatable :: tmp_int(:)
real(kind=REAL64),allocatable :: fr(:)
integer :: NumCol
integer, allocatable :: iTable_(:,:)
real(kind=REAL64), allocatable :: rTable_(:,:)
call formatter%open(File, pFIO_READ, rc=status)
meta = formatter%read(rc=status)
ref => meta%get_attribute('N_Grids')
attr_val => ref%get_values()
select type(attr_val)
type is (integer(INT32))
ng = attr_val(1)
endselect
if (present(n_Grids)) then
n_Grids = ng
endif
ntile = meta%get_dimension('tile')
if (present(n_tiles)) then
n_tiles = ntile
endif
if (present(nx)) then
ref => meta%get_attribute('raster_nx')
attr_val => ref%get_values()
select type(attr_val)
type is (integer(INT32))
nx = attr_val(1)
endselect
endif
if (present(ny)) then
ref => meta%get_attribute('raster_ny')
attr_val => ref%get_values()
select type (attr_val)
type is (integer(INT32))
ny = attr_val(1)
endselect
endif
if (present(N_PfafCat)) then
ref => meta%get_attribute('N_PfafCat')
attr_val => ref%get_values()
select type (attr_val)
type is (integer(INT32))
N_PfafCat = attr_val(1)
endselect
endif
do ll = 1, ng
if (ll == 1) then
ocn_str = ''
else
ocn_str = '_ocn'
endif
if (present(GridName)) then
attr = 'Grid'//trim(ocn_str)//'_Name'
ref =>meta%get_attribute(attr)
char_val => ref%get_value()
select type(char_val)
type is(character(*))
GridName(ll) = char_val
class default
print('unsupported subclass (not string) of attribute named '//attr)
end select
endif
if (present(IM)) then
attr = 'IM'//trim(ocn_str)
ref =>meta%get_attribute(attr)
attr_val => ref%get_values()
select type(attr_val)
type is( integer(INT32))
IM(ll) = attr_val(1)
end select
endif
if (present(JM)) then
attr = 'JM'//trim(ocn_str)
ref =>meta%get_attribute(attr)
attr_val => ref%get_values()
select type(attr_val)
type is(integer(INT32))
JM(ll) = attr_val(1)
end select
endif
enddo
if (present(iTable) .or. present(AVR) ) then
allocate(iTable_(ntile,0:7))
allocate(tmp_int(ntile))
call formatter%get_var('typ', iTable_(:,0))
do ll = 1, ng
if (ll == 1) then
ocn_str = ''
else
ocn_str = '_ocn'
endif
call formatter%get_var('i_indg' //trim(ocn_str), tmp_int, rc=status)
iTable_(:,ll*2) = tmp_int
call formatter%get_var('j_indg' //trim(ocn_str), tmp_int, rc=status)
iTable_(:,ll*2+1) = tmp_int
call formatter%get_var('dummy_index'//trim(ocn_str), tmp_int, rc=status)
if ( ng == 1) then
iTable_(:,4) = tmp_int ! for ease, it is pfaf
! set this 7th column to 1. This is to reproduce a potential bug
! when it is ease grid and mask file is not GEOS5_10arcsec_mask
iTable_(:,7) = 1
else
iTable_(:,5+ll) = tmp_int
endif
enddo
call formatter%get_var('pfaf_index', tmp_int, rc=status)
if (ng == 2) then
where (iTable_(:,0) == 100)
iTable_(:,4) = tmp_int
endwhere
endif
endif
if (present(rTable) .or. present(AVR) ) then
allocate(rTable_(ntile,10))
call formatter%get_var('com_lon', rTable_(:,1), rc=status)
call formatter%get_var('com_lat', rTable_(:,2), rc=status)
call formatter%get_var('area', rTable_(:,3), rc=status)
do ll = 1, ng
if (ll == 1) then
ocn_str = ''
else
ocn_str = '_ocn'
endif
call formatter%get_var('frac_cell' //trim(ocn_str), rTable_(:,3+ll), rc=status)
enddo
call formatter%get_var('min_lon', rTable_(:, 6), rc=status)
call formatter%get_var('max_lon', rTable_(:, 7), rc=status)
call formatter%get_var('min_lat', rTable_(:, 8), rc=status)
call formatter%get_var('max_lat', rTable_(:, 9), rc=status)
call formatter%get_var('elev', rTable_(:,10), rc=status)
endif
if (present(AVR)) then
! In GEOSgcm, it already assumes ng = 2, so NumCol = 10
NumCol = NumGlobalVars+NumLocalVars*ng
allocate(AVR(ntile, NumCol))
AVR(:, 1) = iTable_(:,0)
! for EASE grid, the second collum is replaced by the area
AVR(:, 2) = rTable_(:,3)
AVR(:, 3) = rTable_(:,1)
AVR(:, 4) = rTable_(:,2)
AVR(:, 5) = iTable_(:,2)
AVR(:, 6) = iTable_(:,3)
AVR(:, 7) = rTable_(:,4)
if (ng == 1) then
AVR(:,8) = iTable_(:,4)
else
AVR(:, 8) = iTable_(:,6)
AVR(:, 9) = iTable_(:,4)
AVR(:, 10) = iTable_(:,5)
AVR(:, 11) = rTable_(:,5)
AVR(:, 12) = iTable_(:,7)
endif
endif
if (present(iTable)) then
call move_alloc(iTable_, iTable)
endif
if (present(rTable)) then
call move_alloc(rTable_, rTable)
do ll = 1, ng
where ( rTable(:,3+ll) /=0.0 ) rTable(:,3+ll) = rTable(:,3)/rTable(:,3+ll)
enddo
endif
_RETURN(_SUCCESS)
end subroutine MAPL_ReadTilingNC4