ESMF_UGridGetCoords Subroutine

subroutine ESMF_UGridGetCoords (filename, meshid, coords,  &
                                start, count, centerflag, rc)

    character(len=*), intent(in)       :: filename
    integer, intent(in)               :: meshid
    real(ESMF_KIND_R8), pointer        :: coords (:,:)
    integer,  intent(in)               :: count
    integer,  intent(in)               :: start
    logical, intent(in)                :: centerflag
    integer, intent(out)               :: rc

    type(ESMF_VM) :: vm
    integer :: PetNo, PetCnt

    integer :: ncid, VarId
    integer :: ncStatus

    integer :: i, coorddims,pos1,pos2,len
    character(len=256) :: units, errmsg, attname, coordString
    character(len=256), pointer :: coordNames(:)
    real(ESMF_KIND_R8) :: earthradius

#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 count==0, return with success
    if (count==0) then
       rc=ESMF_SUCCESS
       return
    endif
    ncStatus = nf90_open (path=trim(filename), mode=nf90_nowrite, ncid=ncid)
    if (CDFCheckError (ncStatus, &
      ESMF_METHOD,  &
      ESMF_SRCLINE, trim(filename), &
      rc)) return

    coorddims = size(coords, 2)
    if (centerflag) then
       attname="face_coordinates"
    else
       attname="node_coordinates"
    endif
    ncStatus = nf90_inquire_attribute(ncid, meshId, attname, len=len)
    errmsg = "Attribute "//trim(attname)//" in "//trim(filename)
    if (CDFCheckError (ncStatus, &
          ESMF_METHOD, &
          ESMF_SRCLINE,&
          errmsg,&
          rc)) return
    ncStatus = nf90_get_att (ncid, meshId, attname, coordString)
    if (CDFCheckError (ncStatus, &
      ESMF_METHOD,  &
      ESMF_SRCLINE, errmsg, &
      rc)) return

    ! Parse the attribute to find the varible names for node_coordinates
    if (coorddims == 2) then
       allocate( coordNames(2))
       pos1 = index(coordString(1:)," ")
       coordNames(1) = trim(coordString(1:pos1-1))
       coordNames(2)=trim(coordString(pos1+1:len))
    else
       allocate( coordNames(3))
       pos1 = index(coordString(1:)," ")
       coordNames(1) = trim(coordString(1:pos1-1))
       pos2 = index(coordString(pos1+1:)," ")
       coordNames(2) = trim(coordString(pos1+1:pos1+pos2-1))
       coordNames(3)=trim(coordString(pos1+pos2+1:len))
    endif

    do i=1,coorddims
      errmsg = "Variable "//trim(coordNames(i))//" in "//trim(filename)
      ncStatus = nf90_inq_varid (ncid, trim(coordNames(i)), VarId)
      if (CDFCheckError (ncStatus, &
        ESMF_METHOD,  &
        ESMF_SRCLINE, errmsg, &
        rc)) return
      ncStatus = nf90_get_var (ncid, VarId, coords(:,i) , start=(/start/), &
                         count=(/count/))
      if (CDFCheckError (ncStatus, &
        ESMF_METHOD,  &
        ESMF_SRCLINE, errmsg, &
        rc)) return

      ! Convert to Cartisian 3D coordinates
      ! get the attribute 'units'
      ncStatus = nf90_inquire_attribute(ncid, VarId, "units", len=len)
      errmsg = "Attribute units for "//coordNames(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
            coords(:,i) = coords(:,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
        coords(:,i)=1+coords(:,i)/earthradius
      endif
    enddo
    ncStatus = nf90_close(ncid)
    if (CDFCheckError (ncStatus, &
        ESMF_METHOD,  &
        ESMF_SRCLINE, trim(filename), &
        rc)) return
#else
    call ESMF_LogSetError(rcToCheck=ESMF_RC_LIB_NOT_PRESENT, &
                 msg="- ESMF_NETCDF not defined when lib was compiled", &
                 ESMF_CONTEXT, rcToReturn=rc)
    return
#endif

end subroutine ESMF_UGridGetCoords