function ESMF_MeshCreateCubedSphere(tileSize, nx, ny, name, rc)
! !RETURN VALUE:
type(ESMF_Mesh) :: ESMF_MeshCreateCubedSphere
! !ARGUMENTS:
integer, intent(in) :: tileSize
integer, intent(in) :: nx
integer, intent(in) :: ny
character(len=*), intent(in), optional :: name
integer, intent(out), optional :: rc
!
! !DESCRIPTION:
! Create a {\tt ESMF\_Mesh} object for a cubed sphere grid using identical regular decomposition for every tile.
! The grid coordinates are generated based on the algorithm used by GEOS-5, The tile resolution is defined by
! {\tt tileSize}. Each tile is decomposed into nx x ny blocks and the total number of DEs used
! is nx x ny x 6. If the total PET is not equal to the number of DEs, the DEs are distributed
! into PETs in the default cyclic distribution. Internally, the nodes and the elements from multiple DEs are
! collapsed into a 1D array. Therefore, the nodal distgrid or the element distgrid attached to the Mesh object
! is always a one DE arbitrarily distributed distgrid. The sequential indices of the nodes and the elements
! are derived based on the location of the point in the Cubed Sphere grid. If an element is located at {\tt (x, y)} of
! tile {\tt n}. Its sequential index would be {\tt (n-1)*tileSize*tileSize+(y-1)*tileSize+x}. If it is a node, its
! sequential index would be {\tt (n-1)*(tileSize+1)*(tileSize+1)+(y-1)*(tileSize+1)+x}.
!
! The arguments are:
! \begin{description}
! \item[tilesize]
! The number of elements on each side of the tile of the Cubed Sphere grid
! \item[nx]
! The number of blocks on the horizontal size of each tile
! \item[ny]
! The number of blocks on the vertical size of each tile
! \item [{[name]}]
! The name of the Mesh.
! \item[{[rc]}]
! Return code; equals {\tt ESMF\_SUCCESS} if there are no errors.
! \end{description}
!
!EOP
type(ESMF_Mesh) :: mesh
type(ESMF_VM) :: vm
integer :: PetNo, PetCnt
integer, parameter :: f_p = selected_real_kind(15) ! double precision
real(ESMF_KIND_R8), allocatable :: lonEdge(:,:), latEdge(:,:)
real(ESMF_KIND_R8), allocatable :: lonEdge1D(:), latEdge1D(:)
real(ESMF_KIND_R8), allocatable :: NodeCoords(:), CenterCoords(:)
real(ESMF_KIND_R8), allocatable :: lonCenter(:,:), latCenter(:,:)
integer, allocatable :: ElemIds(:), NodeIds(:)
integer, allocatable :: ElemType(:)
integer, allocatable :: ElemConn(:)
integer, allocatable :: NodeOwners(:)
integer :: sizei, sizej, starti, startj, tile
integer :: rem, rem1, rem2, ind
integer :: i, j, k, kk, kksave, l
integer :: localNodes, localElems
integer :: totalNodes
integer, allocatable :: firstOwners(:), recvbuf(:), map(:)
integer, allocatable :: origIds(:)
integer :: maxDuplicate, uniquenodes
integer :: localrc
real(ESMF_KIND_R8) :: start_lat, end_lat, TOL
real(ESMF_KIND_R8) :: starttime, endtime
character(len=80) :: filename1
integer, allocatable :: start(:,:), count(:,:)
integer, allocatable :: DElist(:), tileno(:)
integer :: totalDE, localDE, unqlocalNodes
integer, allocatable :: GlobalIDs(:), localIDs(:)
! 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
#if 0
if (nx * ny * 6 /= PetCnt) then
call ESMF_LogSetError(ESMF_RC_ARG_WRONG, &
msg="nx * ny does not equal to the total number of PETs", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
#endif
allocate(lonEdge(tileSize+1, (tileSize+1)*6),latEdge(tileSize+1, (tileSize+1)*6))
! Create a mesh
mesh = ESMF_MeshCreate(2, 2, coordSys=ESMF_COORDSYS_SPH_DEG, name=name, rc=localrc)
! Distribute center coordinates according to the nx/ny decomposition
#if 0
xsize = tileSize/nx
ysize = tileSize/ny
starty = ysize*(PetNo/nx)+1
startx = xsize*mod(PetNo,nx)+1
tile = (starty-1)/tilesize
#endif
! use actual cubed sphere coordiantes
totalDE = nx*ny*6
localNodes = 0
localElems = 0
if (PetCnt >= totalDE) then
if (PetNo < totalDE) then
localDE=1
else
localDE=0
endif
else
! multiple DEs per PET
localDE=totalDE/PetCnt
rem = mod(totalDE, PetCnt)
if (PetNo < rem) localDE=localDE+1
endif
if (localDE > 0) then
allocate(DElist(localDE), tileno(localDE))
allocate(start(2,localDE), count(2,localDE))
! deNo is zero-based
do i=1,localDE
DElist(i)=PetNo+PetCnt*(i-1)
enddo
do i=1,localDE
tileno(i) = DElist(i)/(nx*ny)+1
rem = mod(DElist(i),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 *, PetNo, DElist(i), 'block:', starti, startj, sizei, sizej, tileno(i)
start(1,i)=starti
start(2,i)=startj
count(1,i)=sizei
count(2,i)=sizej
localNodes = localNodes + (sizei+1)*(sizej+1)
localElems = localElems + sizei*sizej
enddo
! Add elements -- allocate arrays
allocate(ElemIds(localElems), ElemConn(localElems*4), ElemType(localElems))
allocate(centerCoords(localElems*2))
ElemType = ESMF_MESHELEMTYPE_QUAD
k=1
do i=1,localDE
!call ESMF_VMWtime(starttime, rc=rc)
! Generate glocal edge coordinates and local center coordinates
allocate(lonCenter(count(1,i), count(2,i)), latCenter(count(1,i), count(2,i)))
if (i==1) then
call ESMF_UtilCreateCSCoords(tileSize, lonEdge=lonEdge, latEdge=latEdge, &
start=start(:,i), count=count(:,i), &
tile=tileno(i), lonCenter=lonCenter, latCenter=latCenter)
else
call ESMF_UtilCreateCSCoords(tileSize, start=start(:,i), count=count(:,i), &
tile=tileno(i), lonCenter=lonCenter, latCenter=latCenter)
endif
!call ESMF_VMWtime(endtime, rc=rc)
lonCenter = lonCenter * ESMF_COORDSYS_RAD2DEG
latCenter = latCenter * ESMF_COORDSYS_RAD2DEG
do j=1,count(2,i)
do l=1,count(1,i)
ElemIds(k) = (tileno(i)-1)*tilesize*tilesize+(j+start(2,i)-2)*(tileSize)+start(1,i)+l-1
centerCoords(k*2-1) = lonCenter(l,j)
centerCoords(k*2) = latCenter(l,j)
k=k+1
enddo
enddo
deallocate(lonCenter, latCenter)
enddo
totalnodes = (tileSize+1)*(tileSize+1)*6
! convert radius to degrees
lonEdge = lonEdge * ESMF_COORDSYS_RAD2DEG
latEdge = latEdge * ESMF_COORDSYS_RAD2DEG
!Find unique set of node coordinates
allocate(map(totalnodes))
TOL=0.0000000001
start_lat=-91.0
end_lat = 91.0
allocate(lonEdge1D(totalnodes), latEdge1D(totalnodes))
lonEdge1D = reshape(lonEdge, (/totalnodes/))
latEdge1D = reshape(latEdge, (/totalnodes/))
call c_ESMC_ClumpPntsLL(totalNodes, lonEdge1D, latEdge1D, &
TOL, map, uniquenodes, &
maxDuplicate, start_lat, end_lat, rc)
deallocate(lonEdge1D, latEdge1D)
! Create a new array to point the new index back to the original index
allocate(origIds(uniquenodes))
origIds(:)=0
do i=1,totalnodes
k=map(i)+1
if (origIds(k)==0) origIds(k)=i
enddo
! Find total unique local nodes in each PET
allocate(firstowners(totalNodes), recvbuf(totalNodes))
! Global ID for the elements and the nodes using its 2D index: (y-1)*tileSize+x
! for nodes shared by multiple PETs, the PET with smaller rank will own the
! nodes.
k=1
firstOwners = PetCnt+1
!nodeOwners = PetNo
allocate(GlobalIds(localNodes), localIds(localNodes))
do i=1,localDE
do j=(tileno(i)-1)*(tileSize+1)+start(2,i),count(2,i)+start(2,i)+(tileno(i)-1)*(tileSize+1)
do l=start(1,i),count(1,i)+start(1,i)
!use the new index to the unique set of nodes
kk = (j-1)*(tileSize+1)+l
GlobalIds(k) = origIds(map(kk)+1)
k=k+1
enddo
enddo
enddo
! Find if there are any duplicate globalIDs in the array
! First sort GlobalIds, return the original index in localIds and number of unique nodes
call sort_int(GlobalIds, localIds, unqlocalnodes)
allocate(NodeIds(unqlocalNodes), nodeCoords(unqlocalNodes*2), nodeOwners(unqlocalNodes))
do i=1,localnodes
firstOwners(GlobalIds(i))=PetNo
enddo
! global minimum of firstOwners to find the owner of the nodes (use the smallest PetNo)
call ESMF_VMAllReduce(vm, firstOwners, recvbuf, totalNodes, &
ESMF_REDUCE_MIN, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i=1,localNodes
k=localIds(i)
NodeIds(k)=GlobalIds(i)
nodeOwners(k)=recvbuf(GlobalIds(i))
enddo
k=1
kksave = 0
do i=1,localDE
do j=(tileno(i)-1)*(tileSize+1)+start(2,i),count(2,i)+start(2,i)+(tileno(i)-1)*(tileSize+1)
do l=start(1,i),count(1,i)+start(1,i)
kk=localIds(k)
if (kk > kksave) then
nodeCoords(kk*2-1)=lonEdge(l, j)
nodeCoords(kk*2)=latEdge(l,j)
endif
k=k+1
kksave = kk
enddo
enddo
enddo
deallocate(firstOwners, recvbuf, map, origIds)
call ESMF_MeshAddNodes(mesh, NodeIds=NodeIds, &
NodeCoords = NodeCoords, NodeOwners = NodeOwners, &
rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!deallocate(NodeIds, NodeCoords, NodeOwners)
k=1 !local element index
kk=1 !local node index
do i=1,localDE
do j=1,count(2,i)
do l=1,count(1,i)
ElemConn(k*4-3)=localIds(kk)
ElemConn(k*4-2)=localIds(kk+1)
ElemConn(k*4-1)=localIds(kk+1+(count(1,i)+1))
ElemConn(k*4)=localIds(kk+(count(1,i)+1))
k=k+1
kk=kk+1
enddo
kk=kk+1
enddo
kk=kk+count(1,i)+1
enddo
call ESMF_MeshAddElements(mesh, ElemIds, ElemType, ElemConn, &
elementCoords=centerCoords, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
deallocate(ElemIds, ElemConn, ElemType, centerCoords)
else !localDE=0
!still have to call ESMF_MeshAddNodes() and ESMF_MeshAddElements() even if there is no DEs
!First participate in ESMF_VMALlReduce()
totalnodes = (tileSize+1)*(tileSize+1)*6
allocate(firstowners(totalnodes), recvbuf(totalnodes))
firstOwners = PetCnt+1
! global minimum of firstOwners to find the owner of the nodes (use the smallest PetNo)
call ESMF_VMAllReduce(vm, firstOwners, recvbuf, totalnodes, &
ESMF_REDUCE_MIN, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
allocate(NodeIds(0), NodeCoords(0), NodeOwners(0))
call ESMF_MeshAddNodes(mesh, NodeIds=NodeIds, &
NodeCoords = NodeCoords, NodeOwners = NodeOwners, &
rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
allocate(ElemIds(0), ElemType(0), ElemConn(0))
allocate(centerCoords(0))
call ESMF_MeshAddElements(mesh, ElemIds, ElemType, ElemConn, &
elementCoords=centerCoords, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
deallocate(firstowners, recvbuf)
endif ! localDE>0
deallocate(lonEdge, latEdge)
ESMF_MeshCreateCubedSphere = mesh
! Set return code
if (present(rc)) rc=ESMF_SUCCESS
end function ESMF_MeshCreateCubedSphere