subroutine ESMF_GetNodeFromUGridFile (filename, meshname, nodeCoords, &
nodeCount, startNode, convertToDeg, rc)
character(len=*), intent(in) :: filename
character(len=*), intent(in) :: meshname
real(ESMF_KIND_R8), pointer :: nodeCoords (:,:)
integer, intent(in), optional :: nodeCount
integer, intent(in), optional :: startNode
logical, intent(in), optional :: convertToDeg
integer, intent(out), optional :: rc
type(ESMF_VM) :: vm
integer :: PetNo, PetCnt
integer :: ncid, meshid
integer :: DimIds(2), VarId
integer :: ncStatus
integer :: coordinateDims(3), coordDim, meshDim
integer :: i, j, count, totalCounts, MaxNodePerElmt, localFillValue
integer :: localCount, localStart
character(len=256) :: errmsg, locations, locNames(3), elmtConnName
character(len=256) :: nodeCoordString
character(len=80), allocatable :: nodeCoordNames(:)
integer :: pos0, pos, pos1, pos2, n, yesNode
integer :: len, indexBase
integer :: totalnodes
character(len=24) :: units, attbuf
real(ESMF_KIND_R8) :: deg2rad, earthradius
real(ESMF_KIND_R8) :: rad2deg
real(ESMF_KIND_R8) :: coord(3)
logical :: convertToDegLocal
integer, parameter :: nf90_noerror = 0
real(ESMF_KIND_R8), allocatable:: nodeCoord1D(:)
#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
convertToDegLocal = .false.
if (present(convertToDeg)) convertToDegLocal = convertToDeg
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
! 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
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
if (meshDim == 2) then
allocate( nodeCoordNames(2))
pos1 = index(nodeCoordString(1:)," ")
nodeCoordNames(1) = nodeCoordString(1:pos1-1)
nodeCoordNames(2)=nodeCoordString(pos1+1:len)
else
allocate( nodeCoordNames(3))
pos1 = index(nodeCoordString(1:)," ")
nodeCoordNames(1) = nodeCoordString(1:pos1-1)
pos2 = index(nodeCoordString(pos1+1:)," ")
nodeCoordNames(2) = nodeCoordString(pos1+1:pos1+pos2-1)
nodeCoordNames(3)=nodeCoordString(pos1+pos2+1:len)
endif
! print *, pos1, pos2, trim(nodeCoordNames(1)), ' ', trim(nodeCoordNames(2)), ' ', trim(nodeCoordNames(3))
! 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=totalNodes)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
if (present(startNode)) then
localStart = startNode
else
localStart = 1
endif
if (present(nodeCount)) then
if (localStart+nodeCount-1 > totalNodes) then
call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, &
msg="- startNode+nodeCount > node dimension", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
localCount=nodeCount
else
localCount=totalNodes
endif
allocate( nodeCoords(meshDim,localCount), nodeCoord1D(localCount) )
do i=1,meshDim
errmsg = "Variable "//trim(nodeCoordNames(i))//" in "//trim(filename)
ncStatus = nf90_inq_varid (ncid, trim(nodeCoordNames(i)), VarId)
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
ncStatus = nf90_get_var (ncid, VarId, nodeCoord1D, start=(/localStart/), &
count=(/localCount/))
if (CDFCheckError (ncStatus, &
ESMF_METHOD, &
ESMF_SRCLINE, errmsg, &
rc)) return
do j=1,localCount
nodeCoords(i,j)=nodeCoord1D(j)
enddo
! Convert to Cartisian 3D coordinates
! get the attribute 'units'
ncStatus = nf90_inquire_attribute(ncid, VarId, "units", len=len)
errmsg = "Attribute units for "//trim(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(len:len) .eq. achar(0)) len = len-1
if (i==1 .or. i==2) then
! if units is not "degrees" or "radians" return errors
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 (meshDim == 2 .and. convertToDegLocal) then
if (units(1:7) .eq. "radians") then
rad2deg = 180.0/3.141592653589793238
!print *, 'Convert radians to degree ', rad2deg
nodeCoords(i,:) = nodeCoords(i,:)*rad2deg
endif
elseif (meshDim == 3) then
! if units is "degrees", convert it to radians
if (units(1:7) .eq. "degrees") then
deg2rad = 3.141592653589793238/180.0
nodeCoords(i,:) = nodeCoords(i,:)*deg2rad
endif
endif
else
! normalize the height using the earth radius
if (units(1:len) .eq. "meters") then
earthradius = 6371000.0
else if (units(1:len) .eq. "km" .or. units(1:len) .eq. "kilometers") then
earthradius = 6371.0
else
call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, &
msg="- units attribute for height is not meters, km, or kilometers", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
nodeCoords(i,:)=1+nodeCoords(i,:)/earthradius
endif
enddo
! keep the coordinates as Spherical Degree
#if 0
if (meshDim == 3) then
! Convert the coordinates into Cartesian 3D
do i=1,localCount
coord(1)=nodeCoords(3,i)*cos(nodeCoords(1,i))*cos(nodeCoords(2,i))
coord(2)=nodeCoords(3,i)*sin(nodeCoords(1,i))*cos(nodeCoords(2,i))
coord(3)=nodeCoords(3,i)*sin(nodeCoords(2,i))
nodeCoords(:,i)=coord(:)
enddo
endif
#endif
! 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_GetNodeFromUGridFile