subroutine ESMF_GetElemFromUGridFile (filename, meshname, elmtConn, &
elmtNums, startElmt, rc)
character(len=*), intent(in) :: filename
character(len=*), intent(in) :: meshname
integer(ESMF_KIND_I4), pointer :: elmtConn (:)
integer(ESMF_KIND_I4), pointer :: elmtNums (:)
integer, intent(out) :: startElmt
integer, intent(out), optional :: rc
type(ESMF_VM) :: vm
integer PetNo, PetCnt
integer :: ncid, meshid
integer(ESMF_KIND_I4), allocatable :: elmtConnT(:,:)
integer :: DimIds(2), VarId
integer :: ncStatus
integer :: coordDim, meshDim
integer :: i, j, count, nodeCount, MaxNodePerElmt, localFillValue
integer :: localCount, remain, elmtCount
character(len=256) :: locations, locNames(3), elmtConnName
integer :: pos0, pos, pos1, pos2, n, yesNode
integer :: len, indexBase
character(len=24) :: units
character(len=256) :: errmsg
character(len=24) :: attbuf
integer :: totalConnections
integer, parameter :: nf90_noerror = 0
#ifdef ESMF_NETCDF
! Get VM information
call ESMF_VMGetCurrent(vm, rc=rc)
if (rc /= ESMF_SUCCESS) return
! set up local pet info
call ESMF_VMGet(vm, localPet=PetNo, petCount=PetCnt, rc=rc)
if (rc /= ESMF_SUCCESS) return
ncStatus = nf90_open (path=trim(filename), mode=nf90_nowrite, ncid=ncid)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, trim(filename), &
rc)) return
! get the dummy variable meshname which contains the topology data as attributes
ncStatus = nf90_inq_varid (ncid, trim(meshname), meshId)
errmsg = "Dummy Variable "//trim(meshname)//" does not exist in "//trim(filename)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, trim(meshname), &
rc)) return
! Check if cf_role attribute is set
ncStatus = nf90_get_att (ncid, meshId, "cf_role", values=attbuf)
if (ncStatus /= nf90_noerror) then
ncStatus = nf90_get_att (ncid, meshId, "standard_name", values=attbuf)
errmsg = "Attribute cf_role in "//trim(filename)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
endif
if (attbuf(1:13) .ne. "mesh_topology") then
call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, &
msg="- cf_role attribute is not mesh_topology", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
! Get mesh dimension
ncStatus = nf90_get_att (ncid, meshId, "topology_dimension", values=meshDim)
if (ncStatus /= nf90_noerror) then
ncStatus = nf90_get_att (ncid, meshId, "dimension", values=meshDim)
errmsg = "Attribute topology_dimension or dimension in "//trim(filename)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
endif
if (meshDim == 2) then
! Get element connectivity, if it does not exist, bail out
errmsg = "Attribute face_node_connectivity in "//trim(filename)
ncStatus = nf90_inquire_attribute(ncid, meshId, "face_node_connectivity", len=len)
ncStatus = nf90_get_att (ncid, meshId, "face_node_connectivity", values=elmtConnName)
else if (meshDim == 3) then
! Get element connectivity, if it does not exist, bail out
errmsg = "Attribute volume_node_connectivity in "//trim(filename)
ncStatus = nf90_inquire_attribute(ncid, meshId, "volume_node_connectivity", len=len)
ncStatus = nf90_get_att (ncid, meshId, "volume_node_connectivity", values=elmtConnName)
else
! error message -- wrong dimension
endif
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
! Get element connectivity (UGRID convention is transposed compared to others)
errmsg = "Variable "//elmtConnName(1:len)//" in "//trim(filename)
ncStatus = nf90_inq_varid (ncid, elmtConnName(1:len), VarId)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
! Get elmt conn fill value (if there are triangles mixed with squares, eg)
! _FillValue is optional. It is not needed if all the elements have the same
! number of corner nodes
errmsg = "Attribute "//elmtConnName(1:len)//" _FillValue in "//trim(filename)
ncStatus = nf90_get_att (ncid, VarId, "_FillValue", values=localFillValue)
if (ncStatus /= nf90_noerror) localFillValue = -1
! Get start_index attribute to find out the index base (0 or 1)
ncStatus = nf90_get_att (ncid, VarId, "start_index", values=indexBase)
! if not defined, default to 0-based
if (ncStatus /= nf90_noerror) indexBase = 0
! Get dimensions of element connectivity (transposed) (for allocation)
errmsg = "Dimensions of "//trim(elmtConnName)//" in "//trim(filename)
ncStatus = nf90_inquire_variable (ncid, VarId, dimids=DimIds)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
ncStatus = nf90_inquire_dimension (ncid, DimIds(2), len=elmtCount)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
ncStatus = nf90_inquire_dimension (ncid, DimIds(1), len=MaxNodePerElmt)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
! Decompose the element array evenly across all PETs
localCount = elmtCount/PetCnt
remain = mod (elmtCount,PetCnt)
startElmt = localCount*PetNo +1
if (PetNo==PetCnt-1) localCount=localCount+remain
! print *, PetNo, 'Before allocating elmtConn', localCount
allocate(elmtConnT(MaxNodePerElmt,localCount) )
allocate( elmtNums(localCount) )
! Get element connectivity... transposed
! print *, PetNo, 'Before nf90_get_var()', startElmt,localCount
ncStatus = nf90_get_var (ncid, VarId, elmtConnT, start=(/1,startElmt/), &
count=(/MaxNodePerElmt,localCount/))
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
! Get the number of nodes for each element
! if indexBase is 0, change it to 1 based index
totalConnections = 0
elmtNums(:)=MaxNodePerElmt
if (indexBase == 0) then
do i=1,localcount
do j=1,MaxNodePerElmt
! change 0-base to 1-base
if (elmtConnT(j,i) /= localFillValue) then
elmtConnT(j,i)=elmtConnT(j,i)+1
else
! find the first FillValue
elmtNums(i) = j-1
totalConnections = totalConnections+elmtNums(i)
exit
endif
enddo
if (elmtNums(i)==MaxNodePerElmt) then
totalConnections = totalConnections+elmtNums(i)
endif
enddo
else
do i=1,localcount
j = MaxNodePerElmt
do while (elmtConnT(j,i) == localFillValue)
j = j - 1
enddo
elmtNums(i) = j
totalConnections = totalConnections+elmtNums(i)
enddo
endif
allocate(elmtConn(totalConnections))
j=1
do i=1,localcount
elmtConn(j:j+elmtNums(i)-1)=elmtConnT(1:elmtNums(i),i)
j = j + elmtNums(i)
enddo
if (totalConnections /= j-1) then
print *, PetNo, 'totalconnection does match ', totalConnections, j-1
endif
deallocate(elmtConnT)
ncStatus = nf90_close (ncid=ncid)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, trim(filename), &
rc)) return
#else
call ESMF_LogSetError(ESMF_RC_LIB_NOT_PRESENT, &
msg="- ESMF_NETCDF not defined when lib was compiled", &
ESMF_CONTEXT, rcToReturn=rc)
return
#endif
end subroutine ESMF_GetElemFromUGridFile