ESMF_EsmfGetElement Subroutine

public subroutine ESMF_EsmfGetElement(filename, elementConn, elmtNums, startElmt, elementMask, elementArea, centerCoords, convertToDeg, rc)

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: filename
integer(kind=ESMF_KIND_I4), pointer :: elementConn(:)
integer(kind=ESMF_KIND_I4), pointer :: elmtNums(:)
integer, intent(out) :: startElmt
integer(kind=ESMF_KIND_I4), optional, pointer :: elementMask(:)
real(kind=ESMF_KIND_R8), optional, pointer :: elementArea(:)
real(kind=ESMF_KIND_R8), optional, pointer :: centerCoords(:,:)
logical, intent(in), optional :: convertToDeg
integer, intent(out), optional :: rc

Source Code

subroutine ESMF_EsmfGetElement (filename, elementConn, &
                                 elmtNums, startElmt, elementMask, &
                                 elementArea, centerCoords, convertToDeg, rc)

    character(len=*), intent(in)   :: filename
    integer(ESMF_KIND_I4), pointer :: elementConn (:)
    integer(ESMF_KIND_I4), pointer :: elmtNums (:)
#if defined (ESMF_NO_INTEGER_1_BYTE)
    ! TODO: Eventually use F2008 kind 'int8'.
    integer(selected_int_kind(1)), allocatable :: elmtNums_i1(:)
#else
    integer(ESMF_KIND_I1), allocatable :: elmtNums_i1(:)
#endif
    integer,           intent(out) :: startElmt
    integer(ESMF_KIND_I4), pointer, optional :: elementMask (:)
    real(ESMF_KIND_R8), pointer, optional :: elementArea (:)
    real(ESMF_KIND_R8), pointer, optional :: centerCoords (:,:)
    logical, intent(in), optional  :: convertToDeg
    integer, intent(out), optional :: rc

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

    integer :: ncid
    integer :: ncStatus, status
    integer :: RecCnt (2)

    integer :: DimId
    integer :: nodeCnt, ElmtCount, MaxNodePerElmt, coordDim
    integer :: localCount, remain

    integer :: VarNo, VarType
    character(len=256)::errmsg
    character(len=80) :: units
    integer :: len
    logical :: convertToDegLocal
    integer, parameter :: nf90_noerror = 0
    integer :: localPolyBreakValue, startIndex, i, j, k
    logical :: PolyBreakFound
    logical :: isRaggedArray
    integer, allocatable :: elementConnLocal(:,:)
    integer :: totalConn, startConn
    integer :: senddata(1)
    integer, allocatable :: recvdata(:)
    integer :: memstat

#ifdef ESMF_NETCDF

    ! Init variables
     convertToDegLocal = .false.
     isRaggedArray = .false.
    if (present(convertToDeg)) convertToDegLocal = convertToDeg

    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_SHARE, ncid=ncid)
    if (CDFCheckError (ncStatus, &
      ESMF_METHOD,  &
      ESMF_SRCLINE, trim(filename), &
      rc)) return

    ! get number of elmts
    ncStatus = nf90_inq_dimid (ncid, "elementCount", DimId)
    errmsg = "Dimension elementCount in "//trim(filename)
    if (CDFCheckError (ncStatus, &
      ESMF_METHOD,  &
      ESMF_SRCLINE, errmsg, &
      rc)) return

    ncStatus = nf90_inquire_dimension (ncid, DimId, len=ElmtCount)
     if (CDFCheckError (ncStatus, &
      ESMF_METHOD,  &
      ESMF_SRCLINE, errmsg, &
      rc)) return

    ! Get max_verts_per_elmt
    ncStatus = nf90_inq_dimid (ncid, "maxNodePElement", DimId)
    errmsg = "Dimension maxNodePElement in "//trim(filename)
    ! if not exist, check if connectionCount dimension exists, if so it is ragged array
    if (ncStatus /= nf90_noerror) then
        ncStatus = nf90_inq_dimid (ncid, "connectionCount", DimId)
        errmsg = "Either maxNodePElement or connectionCount does not exist"
        if (CDFCheckError (ncStatus, &
           ESMF_METHOD,  &
           ESMF_SRCLINE, errmsg, &
           rc)) return
        isRaggedArray = .true.
    else
      ncStatus = nf90_inquire_dimension (ncid, DimId, len=MaxNodePerElmt)
      if (CDFCheckError (ncStatus, &
        ESMF_METHOD,  &
        ESMF_SRCLINE, errmsg, &
        rc)) return
    endif

    ! Decompose the elmt array evenly on all the PEs
    localcount = elmtCount/PetCnt
    remain = mod(elmtCount,PetCnt)
    startElmt = localcount * PetNo+1
    if (PetNo == (PetCnt-1)) localcount = localcount+remain

    ! allocate memory for elmts
    allocate (elmtNums (localcount), stat=memstat)
    if (ESMF_LogFoundAllocError(memstat,  &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! read num_elmt_verts
    ncStatus = nf90_inq_varid (ncid, "numElementConn", VarNo)
    errmsg = "Variable numElementConn in "//trim(filename)
    if (CDFCheckError (ncStatus, &
      ESMF_METHOD,  &
      ESMF_SRCLINE, errmsg, &
      rc)) return

    ncStatus = nf90_inquire_variable (ncid, VarNo, xtype=VarType)
    errmsg = "Variable numElementConn type inquiry in "//trim(filename)
    if (CDFCheckError (ncStatus, &
      ESMF_METHOD,  &
      ESMF_SRCLINE, errmsg, &
      rc)) return

    ! Even though NetCDF would automatically do data conversion, special casing NF90_BYTE
    ! can mitigate memory issues when reading very large arrays.  In particular, the
    ! Intel compiler can place large temporaries on the stack, rather than heap, causing
    ! problems.  (See ticket 3614272.)
    select case (VarType)
    case (NF90_INT)
      ncStatus = nf90_get_var (ncid, VarNo, elmtNums, start=(/startElmt/), count=(/localcount/))
      errmsg = "Reading numElementConn from int variable in " // trim (filename)
      if (CDFCheckError (ncStatus, &
          ESMF_METHOD,  &
          ESMF_SRCLINE, errmsg, &
          rc)) return

    case (NF90_BYTE)
      allocate (elmtNums_i1(localcount), stat=memstat)
      if (ESMF_LogFoundAllocError(memstat,  &
          ESMF_CONTEXT, rcToReturn=rc)) return

      ncStatus = nf90_get_var (ncid, VarNo, elmtNums_i1, start=(/startElmt/), count=(/localcount/))
      errmsg = "Reading numElementConn from byte variable in " // trim (filename)
      if (CDFCheckError (ncStatus, &
          ESMF_METHOD,  &
          ESMF_SRCLINE, errmsg, &
          rc)) return

      elmtNums = elmtNums_i1
      deallocate (elmtNums_i1, stat=memstat)
      if (ESMF_LogFoundDeallocError(memstat,  &
          ESMF_CONTEXT, rcToReturn=rc)) return

    case default
      if (ESMF_LogFoundError(ESMF_RC_FILE_UNEXPECTED,  &
          msg='unsupport numElementConn variable type', &
          ESMF_CONTEXT, rcToReturn=rc)) return
    end select

    ! calculate the RaggedArray size
    totalConn = 0
    do i=1,localcount
      totalConn=totalConn+elmtNums(i)
    enddo
    allocate(elementConn(totalConn), stat=memstat)
    if (ESMF_LogFoundAllocError(memstat,  &
        ESMF_CONTEXT, rcToReturn=rc)) return

    ! read elmt_verts data
    ncStatus = nf90_inq_varid (ncid, "elementConn", VarNo)
    errmsg = "Variable elementConn in "//trim(filename)
    if (CDFCheckError (ncStatus, &
       ESMF_METHOD,  &
       ESMF_SRCLINE, errmsg, &
       rc)) return

    if (isRaggedArray) then
      ! need to find out the start index, broadcast the totalConn
      senddata(1)=totalConn
      allocate(recvdata(PetCnt), stat=memstat)
      if (ESMF_LogFoundAllocError(memstat,  &
          ESMF_CONTEXT, rcToReturn=rc)) return

      call ESMF_VMAllGather(vm, senddata, recvdata, 1, rc=status)
      if (ESMF_LogFoundError(status, ESMF_ERR_PASSTHRU, &
                ESMF_CONTEXT, rcToReturn=rc)) return
      startConn=1
      do i=1,PetNo
        startConn=startConn+recvdata(i)
      enddo
      deallocate(recvdata, stat=memstat)
      if (ESMF_LogFoundDeallocError(memstat,  &
          ESMF_CONTEXT, rcToReturn=rc)) return

      ncStatus = nf90_get_var(ncid, VarNo, elementConn, start=(/startConn/), &
                              count=(/totalConn/))
      if (CDFCheckError (ncStatus, &
          ESMF_METHOD,  &
          ESMF_SRCLINE, errmsg, &
          rc)) return
    else
       ! allocate local 2D array to read the data
       allocate (elementConnLocal (MaxNodePerElmt, localcount), stat=memstat)
       if (ESMF_LogFoundAllocError(memstat,  &
           ESMF_CONTEXT, rcToReturn=rc)) return

       ncStatus = nf90_get_var(ncid, VarNo, elementConnLocal, start=(/1,startElmt/), &
                              count=(/MaxNodePerElmt, localcount/))
       if (CDFCheckError (ncStatus, &
          ESMF_METHOD,  &
          ESMF_SRCLINE, errmsg, &
          rc)) return
       ! copy it to elementConn
       j=1
       do i=1,localcount
          elementConn(j:j+elmtNums(i)-1)=elementConnLocal(1:elmtNums(i),i)
          j = j + elmtNums(i)
       end do
       deallocate(elementConnLocal, stat=memstat)
       if (ESMF_LogFoundAllocError(memstat,  &
           ESMF_CONTEXT, rcToReturn=rc)) return

    endif

    ! Get polybreak if it exists
    PolyBreakFound=.false.
     ncStatus = nf90_get_att (ncid, VarNo, "polygon_break_value", &
         values=localPolyBreakValue)
    if (ncStatus == nf90_noerror) then
       PolyBreakFound=.true.
    endif

    ! If a Polygon break value was found, then change those to internal polygon break
    if (PolyBreakFound) then
       do i=1,totalConn
          if (elementConn(i)==localPolyBreakValue) then
             elementConn(i)=ESMF_MESH_POLYBREAK
          endif
       enddo
    endif

    ! Get start_index attribute if it exists
    ! change it to 1-based if it is zero-based
    ncStatus = nf90_get_att (ncid, VarNo, "start_index", &
         values=startIndex)
    if (ncStatus == nf90_noerror) then
      if (startIndex == 0) then
         do i=1,totalConn
           if (elementConn(i) >= 0) then
              elementConn(i)=elementConn(i)+1
           endif
         enddo
      endif
    endif

    ! Check for negative index values that is not defiend as ESMF_MESH_POLYBREAK
    do i=1,totalConn
       if (elementConn(i) < 0 .and. elementConn(i) /= ESMF_MESH_POLYBREAK) then
           call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, &
                msg="- negative index found in elementConn table", &
                ESMF_CONTEXT, rcToReturn=rc)
           return
       endif
    enddo          
    if (present(elementMask)) then
       allocate(elementMask(localcount), stat=memstat)
       if (ESMF_LogFoundAllocError(memstat,  &
           ESMF_CONTEXT, rcToReturn=rc)) return

       ncStatus = nf90_inq_varid (ncid, "elementMask", VarNo)
       errmsg = "Variable elementMask in "//trim(filename)
       if (CDFCheckError (ncStatus, &
          ESMF_METHOD,  &
          ESMF_SRCLINE, errmsg, &
          rc)) return

       ncStatus = nf90_get_var (ncid, VarNo, elementMask, start=(/startElmt/), count=(/localcount/))
       if (CDFCheckError (ncStatus, &
          ESMF_METHOD,  &
           ESMF_SRCLINE, errmsg, &
          rc)) return

    end if

    if (present(elementArea)) then
       allocate(elementArea(localcount), stat=memstat)
       if (ESMF_LogFoundAllocError(memstat,  &
           ESMF_CONTEXT, rcToReturn=rc)) return

       ncStatus = nf90_inq_varid (ncid, "elementArea", VarNo)
       errmsg = "Variable elementArea in "//trim(filename)
       if (CDFCheckError (ncStatus, &
          ESMF_METHOD,  &
          ESMF_SRCLINE, errmsg, &
          rc)) return

       ncStatus = nf90_get_var (ncid, VarNo, elementArea, start=(/startElmt/), &
                                count=(/localcount/))
       if (CDFCheckError (ncStatus, &
          ESMF_METHOD,  &
          ESMF_SRCLINE, errmsg, &
          rc)) return
    end if

    ! If centerCoords exists in file, return information
    if (present(centerCoords)) then
       ncStatus = nf90_inq_varid (ncid, "centerCoords", VarNo)
       if (ncStatus == nf90_noerror) then
          ! Get vertex dimension
          ncStatus = nf90_inq_dimid (ncid, "coordDim", DimId)
          errmsg = "Dimension coordDim in "//trim(filename)
          if (CDFCheckError (ncStatus, &
               ESMF_METHOD,  &
               ESMF_SRCLINE, errmsg, &
               rc)) return

          ncStatus = nf90_inquire_dimension (ncid, DimId, len=coordDim)
          if (CDFCheckError (ncStatus, &
               ESMF_METHOD,  &
               ESMF_SRCLINE, errmsg, &
               rc)) return

          allocate(centerCoords(coordDim, localcount), stat=memstat)
          if (ESMF_LogFoundAllocError(memstat,  &
               ESMF_CONTEXT, rcToReturn=rc)) return

          ncStatus = nf90_inq_varid (ncid, "centerCoords", VarNo)
          errmsg = "Variable centerCoords in "//trim(filename)
          if (CDFCheckError (ncStatus, &
               ESMF_METHOD,  &
               ESMF_SRCLINE, errmsg, &
               rc)) return

          ncStatus = nf90_get_var (ncid, VarNo, centerCoords, start=(/1,startElmt/), &
               count=(/coordDim, localcount/))
          if (CDFCheckError (ncStatus, &
               ESMF_METHOD,  &
               ESMF_SRCLINE, errmsg, &
               rc)) return

          ! Check units, but only if 2D
          if (coordDim==2) then
             ncStatus = nf90_inquire_attribute(ncid, VarNo, "units", len=len)
             if (CDFCheckError (ncStatus, &
                  ESMF_METHOD, &
                  ESMF_SRCLINE,errmsg,&
                  rc)) return

             ncStatus = nf90_get_att(ncid, VarNo, "units", units)
             if (CDFCheckError (ncStatus, &
                  ESMF_METHOD, &
                  ESMF_SRCLINE,errmsg,&
                  rc)) return
             ! if len != 7, something is wrong, check the value.  If it starts
             ! with Degres/degrees/Radians/radians, ignore the garbage after the
             ! word.  Otherwise, return the whole thing
             if (units(len:len) .eq. achar(0)) len = len-1
             units = ESMF_UtilStringLowerCase(units(1:len))
             if (units(1:len) .ne. 'degrees' .and. &
                  units(1:len) .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 degree
             if (convertToDegLocal) then
                if (units(1:len) .eq. "radians") then
                   centerCoords(:,:) = &
                        centerCoords(:,:)*ESMF_COORDSYS_RAD2DEG
                endif
             endif
          endif
       endif
    end if

    ncStatus = nf90_close (ncid=ncid)
    if (CDFCheckError (ncStatus, &
      ESMF_METHOD,  &
      ESMF_SRCLINE, trim(filename), &
      rc)) return
#else
    if (ESMF_LogFoundError(ESMF_RC_LIB_NOT_PRESENT, &
                msg="- ESMF_NETCDF not defined when lib was compiled", &
                ESMF_CONTEXT, rcToReturn=rc)) return
#endif

end subroutine ESMF_EsmfGetElement