function ESMF_GridCreateCubedSphereReg(tileSize,keywordEnforcer, &
regDecompPTile, decompflagPTile, &
coordSys, coordTypeKind, &
deLabelList, staggerLocList, &
delayout, indexflag, name, transformArgs, rc)
!
! !RETURN VALUE:
type(ESMF_Grid) :: ESMF_GridCreateCubedSphereReg
!
! !ARGUMENTS:
integer, intent(in) :: tilesize
type(ESMF_KeywordEnforcer), optional:: keywordEnforcer ! must use keywords below
integer, intent(in), optional :: regDecompPTile(:,:)
type(ESMF_Decomp_Flag), target, intent(in), optional :: decompflagPTile(:,:)
type(ESMF_CoordSys_Flag), intent(in), optional :: coordSys
type(ESMF_TypeKind_Flag), intent(in), optional :: coordTypeKind
integer, intent(in), optional :: deLabelList(:)
type(ESMF_StaggerLoc), intent(in), optional :: staggerLocList(:)
type(ESMF_DELayout), intent(in), optional :: delayout
type(ESMF_Index_Flag), intent(in), optional :: indexflag
character(len=*), intent(in), optional :: name
type(ESMF_CubedSphereTransform_Args), intent(in), optional :: transformArgs
integer, intent(out), optional :: rc
!
! !DESCRIPTION:
! Create a six-tile {\tt ESMF\_Grid} for a Cubed Sphere grid using regular decomposition. Each tile can
! have different decomposition. The grid coordinates are generated based on the algorithm used by GEOS-5,
! The tile resolution is defined by tileSize.
!
! The arguments are:
! \begin{description}
! \item[tilesize]
! The number of elements on each side of the tile of the Cubed Sphere grid.
! \item[{[regDecompPTile]}]
! List of DE counts for each dimension. The second index steps through
! the tiles. The total {\tt deCount} is determined as the sum over
! the products of {\tt regDecompPTile} elements for each tile.
! By default every tile is decomposed in the same way. If the total
! PET count is less than 6, one tile will be assigned to one DE and the DEs
! will be assigned to PETs sequentially, therefore, some PETs may have
! more than one DE. If the total PET count is greater than 6, the total
! number of DEs will be a multiple of 6 and less than or equal to the total
! PET count. For instance, if the total PET count is 16, the total DE count
! will be 12 with each tile decomposed into 1x2 blocks. The 12 DEs are mapped
! to the first 12 PETs and the remaining 4 PETs have no DEs locally, unless
! an optional {\tt delayout} is provided.
! \item[{[decompflagPTile]}]
! List of decomposition flags indicating how each dimension of each
! tile is to be divided between the DEs. The default setting
! is {\tt ESMF\_DECOMP\_BALANCED} in all dimensions for all tiles.
! See section \ref{const:decompflag} for a list of valid decomposition
! flag options. The second index indicates the tile number.
! \item[{[deLabelList]}]
! List assigning DE labels to the default sequence of DEs. The default
! sequence is given by the column major order of the {\tt regDecompPTile}
! elements in the sequence as they appear following the tile index.
! \item[{[staggerLocList]}]
! The list of stagger locations to fill with coordinates. Only {\tt ESMF\_STAGGERLOC\_CENTER}
! and {\tt ESMF\_STAGGERLOC\_CORNER} are supported. If not present, no coordinates
! will be added or filled.
! \item[{[coordSys]}]
! The coordinate system of the grid coordinate data.
! Only ESMF\_COORDSYS\_SPH\_DEG and ESMF\_COORDSYS\_SPH\_RAD are supported.
! If not specified then defaults to ESMF\_COORDSYS\_SPH\_DEG.
! \item[{[coordTypeKind]}]
! The type/kind of the grid coordinate data. Only ESMF\_TYPEKIND\_R4
! and ESMF\_TYPEKIND\_R8 are supported.
! If not specified then defaults to ESMF\_TYPEKIND\_R8.
! \item[{[delayout]}]
! Optional {\tt ESMF\_DELayout} object to be used. By default a new
! DELayout object will be created with as many DEs as there are PETs,
! or tiles, which ever is greater. If a DELayout object is specified,
! the number of DEs must match {\tt regDecompPTile}, if present. In the
! case that {\tt regDecompPTile} was not specified, the {\tt deCount}
! must be at least that of the default DELayout. The
! {\tt regDecompPTile} will be constructed accordingly.
! \item[{[indexflag]}]
! Indicates the indexing scheme to be used in the new Grid. Please see
! Section~\ref{const:indexflag} for the list of options. If not present,
! defaults to ESMF\_INDEX\_DELOCAL.
! \item[{[name]}]
! {\tt ESMF\_Grid} name.
! \item[{[transformArgs]}]
! A data type containing the stretch factor, target longitude, and target latitude
! to perform a Schmidt transformation on the Cubed-Sphere grid. See section
! \ref{sec:usage:cubedspherewttransform} for details.
! \item[{[rc]}]
! Return code; equals {\tt ESMF\_SUCCESS} if there are no errors.
! \end{description}
!
!EOP
type(ESMF_VM) :: vm
integer :: PetNo, PetCnt
type(ESMF_DELayout) :: defaultDELayout
type(ESMF_Grid) :: grid, newgrid
type(ESMF_DistGrid) :: distgrid, newdistgrid
integer :: localrc
integer :: i, j
integer :: nx, ny, nxy, bigFac, totalDE
integer :: localDeCount, localDe, DeNo, tile
integer, pointer :: minIndexPTile(:,:)
integer, pointer :: maxIndexPTile(:,:)
type(ESMF_DistGridConnection), pointer :: connectionList(:)
real(kind=ESMF_KIND_R8), pointer :: lonPtrR8(:,:), latPtrR8(:,:)
real(kind=ESMF_KIND_R8), pointer :: lonCornerPtrR8(:,:), latCornerPtrR8(:,:)
real(kind=ESMF_KIND_R4), pointer :: lonPtrR4(:,:), latPtrR4(:,:)
real(kind=ESMF_KIND_R4), pointer :: lonCornerPtrR4(:,:), latCornerPtrR4(:,:)
integer :: tileCount
integer :: starti, startj, sizei, sizej
integer :: ind, rem, rem1, rem2
integer :: start(2), count(2)
integer :: shapLon(2), shapLat(2)
integer, allocatable :: regDecomp(:,:)
integer, allocatable :: demap(:)
integer :: decount
type(ESMF_Index_Flag) :: localIndexFlag
type(ESMF_CoordSys_Flag) :: coordSysLocal
type(ESMF_TypeKind_Flag) :: coordTypeKindLocal
integer :: s
logical :: docenter, docorner
!real(ESMF_KIND_R8) :: starttime, endtime
if (present(rc)) rc=ESMF_SUCCESS
!------------------------------------------------------------------------
! get global vm information
!
call ESMF_VMGetCurrent(vm, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! set up local pet info
call ESMF_VMGet(vm, localPet=PetNo, petCount=PetCnt, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!------------------------------------------------------------------------
! default decomposition. The number of DEs has to be multiple of 6.
! If the total PET count is less than 6, some PETs will get more than one DE.
! Otherwise, total DEs is always less than or equal to total PETs.
#if 1
if (PetCnt < 6) then
totalDE=6
else
totalDE = (PetCnt/6)*6
endif
nxy = totalDE/6
bigFac = 1
do i=2, int(sqrt(float(nxy)))
if ((nxy/i)*i == nxy) then
bigFac = i
endif
enddo
nx = bigFac
ny = nxy/nx
#else
nxy = (PetCnt + 5)/6
totalDE = 6 * nxy
nx = 1
do i = 2, int(sqrt(real(nxy)))
if (mod(nx,i) == 0) nx = i
end do
ny = nxy / nx
#endif
defaultDELayout = ESMF_DELayoutCreate(deCount = totalDE, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
if (present(indexflag)) then
localIndexFlag = indexflag
else
localIndexFlag = ESMF_INDEX_DELOCAL
endif
! Set Default coordSys
if (present(coordSys)) then
if (coordSys .eq. ESMF_COORDSYS_CART) then
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, &
msg="- only ESMF_TYPEKIND_CART is not supported", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
coordSysLocal=coordSys
else
coordSysLocal=ESMF_COORDSYS_SPH_DEG
endif
! Set Default coordTypeKind
if (present(coordTypeKind)) then
if (coordTypeKind .ne. ESMF_TYPEKIND_R4 .and. &
coordTypeKind .ne. ESMF_TYPEKIND_R8) then
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, &
msg="- only ESMF_TYPEKIND_R4 and ESMF_TYPEKIND_R8 are allowed", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
coordTypeKindLocal=coordTypeKind
else
coordTypeKindLocal=ESMF_TYPEKIND_R8
endif
! set defaults
docenter = .false.
docorner = .false.
tileCount = 6
allocate(minIndexPTile(2,tileCount))
allocate(maxIndexPTile(2,tileCount))
allocate(connectionList(12))
call CalculateConnection(tilesize, minIndexPTile, maxIndexPTile, connectionList, localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
allocate(regDecomp(2,6))
regDecomp(1,:)=nx
regDecomp(2,:)=ny
!-------------------------------------------
! - create DistGrid with default decomposition
! must create with ESMF_INDEX_GLOBAL because of how connections were defined
distgrid = ESMF_DistGridCreate(&
minIndexPTile=minIndexPTile, maxIndexPTile=maxIndexPTile, &
regDecompPTile=regDecomp, &
indexflag=ESMF_INDEX_GLOBAL, connectionList=connectionList, &
delayout = defaultDelayout, &
rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
if (present(name)) then
call ESMF_DistGridSet(distgrid, name="DG-"//trim(name), rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
endif
! - create Grid
grid = ESMF_GridCreate(distgrid, coordSys=coordSysLocal, &
coordTypeKind=coordTypeKindLocal, &
indexflag=localIndexFlag, name=name, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
if (present(staggerLocList)) then
call ESMF_DELayoutGet(defaultDElayout, localDeCount = decount, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
if (decount > 0) then
allocate(demap(0:decount-1))
call ESMF_DELayoutGet(defaultDElayout, localDeToDeMap = demap, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
endif
do s=1, size(staggerLocList)
if (staggerLocList(s) == ESMF_STAGGERLOC_EDGE1 .or. &
staggerLocList(s) == ESMF_STAGGERLOC_EDGE2) then
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, &
msg="- only ESMF_STAGGERLOC_CENTER and ESMF_STAGGERLOC_CORNER are supported", &
ESMF_CONTEXT, rcToReturn=rc)
return
elseif (staggerLocList(s) == ESMF_STAGGERLOC_CENTER) then
docenter = .TRUE.
elseif (staggerLocList(s) == ESMF_STAGGERLOC_CORNER) then
docorner = .TRUE.
endif
call ESMF_GridAddCoord(grid, staggerloc=staggerLocList(s), rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
enddo
! calculate the actual cubed sphere coordiantes for each DE
do localDe = 0,decount-1
DeNo = demap(localDe)
tile = DeNo/(nx*ny)+1
rem = mod(DeNo,nx*ny)
sizei = tileSize/nx
sizej = tileSize/ny
rem1 = mod(tileSize, nx)
rem2 = mod(tileSize, ny)
ind = mod(rem,nx)
if (rem1 > 0) then
if (ind < rem1) then
sizei=sizei+1
starti=sizei*ind+1
else
starti=sizei*ind+rem1+1
endif
else
starti = sizei*ind+1
endif
ind = rem/nx
if (rem2 > 0) then
if (ind < rem2) then
sizej=sizej+1
startj=sizej*ind+1
else
startj=sizej*ind+rem2+1
endif
else
startj = sizej*ind+1
endif
!print *, DeNo, 'block:', starti, startj, sizei, sizej, tile
start(1)=starti
start(2)=startj
count(1)=sizei
count(2)=sizej
if (coordTypeKindLocal .eq. ESMF_TYPEKIND_R8) then
if (docenter) then
call ESMF_GridGetCoord(grid, coordDim=1, localDe=localDe, &
farrayPtr=lonPtrR8, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridGetCoord(grid, coordDim=2, localDe=localDe, &
farrayPtr=latPtrR8, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
endif
if (docorner) then
call ESMF_GridGetCoord(grid, coordDim=1, localDe=localDe, &
staggerloc=ESMF_STAGGERLOC_CORNER, farrayPtr=lonCornerPtrR8, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridGetCoord(grid, coordDim=2, localDe=localDe, &
staggerloc=ESMF_STAGGERLOC_CORNER, farrayPtr=latCornerPtrR8, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
endif
else ! ESMF_TYPEKIND_R4
if (docenter) then
call ESMF_GridGetCoord(grid, coordDim=1, localDe=localDe, &
farrayPtr=lonPtrR4, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridGetCoord(grid, coordDim=2, localDe=localDe, &
farrayPtr=latPtrR4, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
allocate(lonPtrR8(count(1), count(2)), latPtrR8(count(1), count(2)))
endif
if (docorner) then
call ESMF_GridGetCoord(grid, coordDim=1, localDe=localDe, &
staggerloc=ESMF_STAGGERLOC_CORNER, farrayPtr=lonCornerPtrR4, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridGetCoord(grid, coordDim=2, localDe=localDe, &
staggerloc=ESMF_STAGGERLOC_CORNER, farrayPtr=latCornerPtrR4, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
shapLon=shape(lonCornerPtrR4) ! make sure lhs and rhs is same shape
shapLat=shape(latCornerPtrR4) ! make sure lhs and rhs is same shape
allocate(lonCornerPtrR8(shapLon(1), shapLon(2)), &
latCornerPtrR8(shapLat(1),shapLat(2)))
endif
endif
!call ESMF_VMWtime(starttime, rc=localrc)
! Generate glocal edge coordinates and local center coordinates
if (docenter .and. docorner) then
call ESMF_UtilCreateCSCoordsPar(tileSize, lonEdge=lonCornerPtrR8, &
latEdge=latCornerPtrR8, start=start, count=count, &
tile=tile, lonCenter=lonPtrR8, latCenter=latPtrR8, &
schmidtTransform=transformArgs)
elseif (docorner) then
call ESMF_UtilCreateCSCoordsPar(tileSize, lonEdge=lonCornerPtrR8, &
latEdge=latCornerPtrR8, start=start, count=count, tile=tile, &
schmidtTransform=transformArgs)
else
call ESMF_UtilCreateCSCoordsPar(tileSize, &
start=start, count=count, &
tile=tile, lonCenter=lonPtrR8, latCenter=latPtrR8, &
schmidtTransform=transformArgs)
endif
!call ESMF_VMWtime(endtime, rc=localrc)
if (coordTypeKindLocal .eq. ESMF_TYPEKIND_R8) then
if (docenter) then
if (coordSysLocal .eq. ESMF_COORDSYS_SPH_DEG) then
lonPtrR8 = lonPtrR8 * ESMF_COORDSYS_RAD2DEG
latPtrR8 = latPtrR8 * ESMF_COORDSYS_RAD2DEG
endif
endif
if (docorner) then
if (coordSysLocal .eq. ESMF_COORDSYS_SPH_DEG) then
lonCornerPtrR8 = lonCornerPtrR8 * ESMF_COORDSYS_RAD2DEG
latCornerPtrR8 = latCornerPtrR8 * ESMF_COORDSYS_RAD2DEG
endif
endif
else ! ESMF_TYPE_KIND_R4
if (docenter) then
if (coordSysLocal .eq. ESMF_COORDSYS_SPH_DEG) then
lonPtrR4 = lonPtrR8 * ESMF_COORDSYS_RAD2DEG
latPtrR4 = latPtrR8 * ESMF_COORDSYS_RAD2DEG
else
lonPtrR4 = lonPtrR8
latPtrR4 = latPtrR8
endif
deallocate(lonPtrR8, latPtrR8)
endif
if (docorner) then
if (coordSysLocal .eq. ESMF_COORDSYS_SPH_DEG) then
lonCornerPtrR4 = lonCornerPtrR8 * ESMF_COORDSYS_RAD2DEG
latCornerPtrR4 = latCornerPtrR8 * ESMF_COORDSYS_RAD2DEG
else
lonCornerPtrR4 = lonCornerPtrR8
latCornerPtrR4 = latCornerPtrR8
endif
deallocate(lonCornerPtrR8, latCornerPtrR8)
endif
endif
!print *, 'Create CS size ', tileSize, 'in', (endtime-starttime)*1000.0, ' msecs'
end do
endif
! Create another distgrid with user specified decomposition
if (present(decompflagPTile) .or. present(regDecompPTile) .or. &
present(delabelList) .or. present(delayout)) then
newdistgrid = ESMF_DistGridCreate(&
minIndexPTile=minIndexPTile, maxIndexPTile=maxIndexPTile, &
regDecompPTile=regDecompPTile, &
decompflagPTile=decompflagPTile, &
delabelList = delabelList, &
indexflag=ESMF_INDEX_GLOBAL, connectionList=connectionList, &
delayout = delayout, &
rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!Redist the grid with the new distgrid
newgrid = ESMF_GridCreate(grid, newdistgrid, name=name, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Destroy old grid
call ESMF_GridDestroy(grid, noGarbage=.true., rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Destroy old distgrid
call ESMF_DistGridDestroy(distgrid, noGarbage=.true., rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
ESMF_GridCreateCubedSphereReg = newgrid
else
ESMF_GridCreateCubedSphereReg = grid
endif
! - deallocate connectionList
deallocate(minIndexPTile, maxIndexPTile)
deallocate(connectionList)
return
end function ESMF_GridCreateCubedSphereReg