subroutine ESMF_GetMesh3DFromUGrid (filename, ncid, meshid, nodeCoords, elmtConn, &
elmtNums, startElmt, faceCoords, 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(:,:)
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
real(ESMF_KIND_R8) :: deg2rad, earthradius
real(ESMF_KIND_R8) :: coord(3)
integer, parameter :: nf90_noerror = 0
real(ESMF_KIND_R8), allocatable:: nodeCoord1D(:), faceCoord1D(:)
integer :: localrc
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
! 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(3))
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
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)
! 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(3) )
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
pos1 = index(faceCoordString(1:)," ")
faceCoordNames(1) = faceCoordString(1:pos1-1)
pos2 = index(faceCoordString(pos1+1:)," ")
faceCoordNames(2)=faceCoordString(pos1+1:pos1+pos2-1)
faceCoordNames(3)=faceCoordString(pos1+pos2+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(3,nodeCount), nodeCoord1D(nodeCount) )
do i=1,3
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
! Convert to Cartisian 3D coordinates
! 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 (i==1 .or. i==2) then
! if units is not "degrees" or "radians" return errors
units = ESMF_UtilStringLowerCase(units(1:len))
if (units(len:len) .eq. achar(0)) len = len-1
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 (units(1:7) .eq. "radians") then
nodeCoords(i,:) = nodeCoords(i,:)*ESMF_COORDSYS_RAD2DEG
endif
else
! normalize the height using the earth radius
if (units(len:len) .eq. achar(0)) len = len-1
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
! save coordinates as ESMF_COORDSYS_SPH_DEG, no need to convert to CART
#if 0
! Convert the coordinates into Cartesian 3D
do i=1,nodeCount
call c_esmc_sphdeg_to_cart(nodeCoords(1,i), nodeCoords(2,i), &
coord(1), coord(2), coord(3), &
localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
nodeCoords(1,i)=nodeCoords(3,i)*coord(1)
nodeCoords(2,i)=nodeCoords(3,i)*coord(2)
nodeCoords(3,i)=nodeCoords(3,i)*coord(3)
enddo
#endif
! 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)
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
! Currently, only hexahedron and tetrahedrons are supported, but
! eventually, we want to support prisms (or wedge as used in UGRID)
errmsg = "Attribute "//elmtConnName(1:len)//" _FillValue in "//trim(filename)
ncStatus = nf90_get_att (ncid, VarId, "_FillValue", values=localFillValue)
if (ncStatus .ne. 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(3,localCount), faceCoord1D(localCount))
do i=1,3
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
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
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
if (totalConnections /= j-1) then
print *, PetNo, 'totalconnection does match ', totalConnections, j-1
endif
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_GetMesh3DFromUGrid