subroutine ESMF_GetMesh2DFromUGrid (filename, ncid, meshid, nodeCoords, elmtConn, &
elmtNums, startElmt, faceCoords, convertToDeg, rc)
character(len=*), intent(in) :: filename
integer, intent(in) :: ncid, meshid
real(ESMF_KIND_R8), pointer :: nodeCoords (:,:)
integer(ESMF_KIND_I4), pointer :: elmtConn (:)
integer(ESMF_KIND_I4), pointer :: elmtNums (:)
integer, intent(out) :: startElmt
real(ESMF_KIND_R8), pointer, optional :: faceCoords(:,:)
logical, intent(in), optional :: convertToDeg
integer, intent(out), optional :: rc
type(ESMF_VM) :: vm
integer PetNo, PetCnt
integer(ESMF_KIND_I4), allocatable :: elmtConnT(:,:)
integer :: DimIds(2), VarId
integer :: ncStatus
integer :: coordinateDims(3), coordDim, meshDim
integer :: i, j, count, nodeCount, MaxNodePerElmt, localFillValue
integer :: localCount, remain, elmtCount, localPolyBreakValue
character(len=256) :: errmsg, locations, locNames(3), elmtConnName
character(len=256) :: nodeCoordString, faceCoordString
character(len=80), allocatable :: nodeCoordNames(:), faceCoordNames(:)
integer :: pos0, pos, pos1, pos2, n, yesNode
integer :: len, indexBase
character(len=24) :: units
logical :: convertToDegLocal
integer, parameter :: nf90_noerror = 0
real(ESMF_KIND_R8), allocatable:: nodeCoord1D(:), faceCoord1D(:)
integer :: totalConnections
#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
if (present(convertToDeg)) then
convertToDegLocal = convertToDeg
else
convertToDegLocal = .false.
endif
! Get node coordinates
ncStatus = nf90_inquire_attribute(ncid, meshId, "node_coordinates", len=len)
errmsg = "Attribute node_coordinates in "//trim(filename)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
allocate( nodeCoordNames(2) )
ncStatus = nf90_get_att (ncid, meshId, "node_coordinates", nodeCoordString)
errmsg = "Attribute node_coordinates in "//trim(filename)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
! Parse the attribute to find the varible names for node_coordinates
pos = index(nodeCoordString(1:)," ")
nodeCoordNames(1) = nodeCoordString(1:pos-1)
nodeCoordNames(2)=nodeCoordString(pos+1:len)
! Get face coordinates if faceCoords argument is given
if (present(faceCoords)) then
ncStatus = nf90_inquire_attribute(ncid, meshId, "face_coordinates", len=len)
errmsg = "Attribute face_coordinates in "//trim(filename)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
allocate( faceCoordNames(2) )
ncStatus = nf90_get_att (ncid, meshId, "face_coordinates", faceCoordString)
errmsg = "Attribute face_coordinates in "//trim(filename)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
! Parse the attribute to find the varible names for face_coordinates
pos = index(faceCoordString(1:)," ")
faceCoordNames(1) = faceCoordString(1:pos-1)
faceCoordNames(2)=faceCoordString(pos+1:len)
endif
! Get dimension (# nodes) used to define the node coordinates
errmsg = "Variable "//trim(nodeCoordNames(1))//" in "//trim(filename)
ncStatus = nf90_inq_varid (ncid, nodeCoordNames(1), VarId)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
ncStatus = nf90_inquire_variable (ncid, VarId, dimids=DimIds)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
ncStatus = nf90_inquire_dimension (ncid, DimIds(1), len=nodeCount)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
allocate( nodeCoords(2,nodeCount), nodeCoord1D(nodeCount) )
do i=1,2
errmsg = "Variable "//nodeCoordNames(i)//" in "//trim(filename)
ncStatus = nf90_inq_varid (ncid, nodeCoordNames(i), VarId)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
ncStatus = nf90_get_var (ncid, VarId, nodeCoord1D)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
do j=1,nodeCount
nodeCoords(i,j)=nodeCoord1d(j)
enddo
! if units is radians_east or radians_north, convert to degrees
! get the attribute 'units'
ncStatus = nf90_inquire_attribute(ncid, VarId, "units", len=len)
errmsg = "Attribute units for "//nodeCoordNames(i)//" in "//trim(filename)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE,&
errmsg,&
rc)) return
ncStatus = nf90_get_att(ncid, VarId, "units", units)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE,&
errmsg,&
rc)) return
! if units is not "degrees" or "radians" return errors
if (units(len:len) .eq. achar(0)) len = len-1
units = ESMF_UtilStringLowerCase(units(1:len))
if (units(1:7) .ne. 'degrees' .and. units(1:7) .ne. 'radians') then
call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, &
msg="- units attribute is not degrees or radians", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
! if units is "radians", convert it to degrees
if (convertToDegLocal) then
if (units(1:7) .eq. "radians") then
nodeCoords(i,:) = &
nodeCoords(i,:)*ESMF_COORDSYS_RAD2DEG
endif
endif
enddo
! 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)
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 PolyBreak Value if it's not present, then set it to localFillValue, so
! it's ignored
ncStatus = nf90_get_att (ncid, VarId, "polygon_break_value", &
values=localPolyBreakValue)
if (ncStatus /= nf90_noerror) then
localPolyBreakValue = localFillValue
else
! If it's been set, then make sure that this value isn't the same as _FillValue
if (localPolyBreakValue == localFillValue) then
call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, &
msg="- polygon_break_value can't be the same as _localFillValue", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
endif
! 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
allocate(elmtConnT(MaxNodePerElmt,localCount) )
allocate( elmtNums(localCount) )
! Get element connectivity... transposed
ncStatus = nf90_get_var (ncid, VarId, elmtConnT, start=(/1,startElmt/), &
count=(/MaxNodePerElmt,localCount/))
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
! Get faceCoordinates if faceCoords is present
if (present(faceCoords)) then
allocate(faceCoords(2,localCount), faceCoord1D(localCount))
do i=1,2
errmsg = "Variable "//faceCoordNames(i)//" in "//trim(filename)
ncStatus = nf90_inq_varid (ncid, faceCoordNames(i), VarId)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
ncStatus = nf90_get_var (ncid, VarId, faceCoord1D, start=(/startElmt/), &
count=(/localCount/))
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
do j=1,localCount
faceCoords(i,j)=faceCoord1d(j)
enddo
! if units is radians_east or radians_north, convert to degrees
! get the attribute 'units'
ncStatus = nf90_inquire_attribute(ncid, VarId, "units", len=len)
errmsg = "Attribute units for "//faceCoordNames(i)//" in "//trim(filename)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE,&
errmsg,&
rc)) return
ncStatus = nf90_get_att(ncid, VarId, "units", units)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE,&
errmsg,&
rc)) return
! if units is not "degrees" or "radians" return errors
if (units(len:len) .eq. achar(0)) len = len-1
units = ESMF_UtilStringLowerCase(units(1:len))
if (units(1:7) .ne. 'degrees' .and. units(1:7) .ne. 'radians') then
call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, &
msg="- units attribute is not degrees or radians", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
! if units is "radians", convert it to degrees
if (convertToDegLocal) then
if (units(1:7) .eq. "radians") then
faceCoords(i,:) = &
faceCoords(i,:)*ESMF_COORDSYS_RAD2DEG
endif
endif
enddo
deallocate(faceCoordNames, faceCoord1d)
endif
! 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
if (elmtConnT(j,i) /= localPolyBreakValue) then
elmtConnT(j,i)=elmtConnT(j,i)+1
endif
else
! find the first FillValue
elmtNums(i) = j-1
totalConnections = totalConnections+elmtNums(i)
exit
endif
enddo
!! Change j to i in elmtNums(i)
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
! Change File PolyBreak to MeshPolyBreak
! (if localPolyBreakValue was set then it'll be different than localFillValue)
if (localPolyBreakValue /= localFillValue) then
do i=1,localcount
do j=1,elmtNums(i)
if (elmtConnT(j,i)==localPolyBreakValue) then
elmtConnT(j,i)=ESMF_MESH_POLYBREAK
endif
enddo
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
deallocate(elmtConnT)
! Deallocations
deallocate( nodeCoordNames, nodeCoord1D )
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_GetMesh2DFromUGrid