readContacts Subroutine

subroutine readContacts(ncid, varid, dims, mosaicname, tilenames, &
                 contact, connindex, rc)

  integer, intent(in)                       :: ncid, varid
  integer, intent(in)                       :: dims(2)
  character(len=*), intent(in)              :: mosaicname
  character(len=*), intent(in), target      :: tilenames(:)
  integer, intent(inout)                     :: contact(:,:)
  integer, intent(inout)                    :: connindex(:,:,:)
  integer, intent(out)                      :: rc

  ! varid is the id for the "contact_region" variable
  ! its standard_name should be "grid_contact_spec" and contact_type should be "boundary"
  ! The "contact_index" attribute defines the variable that defines the edge connection
  ! The contact_region array contacts the neighboring pair of tiles for each connection:
  ! contacts =
  !  "C48_mosaic:tile1::C48_mosaic:tile2",
  !  "C48_mosaic:tile1::C48_mosaic:tile3",
  !  "C48_mosaic:tile1::C48_mosaic:tile5",
  !  "C48_mosaic:tile1::C48_mosaic:tile6",
  !  "C48_mosaic:tile2::C48_mosaic:tile3",
  !  "C48_mosaic:tile2::C48_mosaic:tile4",
  !  "C48_mosaic:tile2::C48_mosaic:tile6",
  !  "C48_mosaic:tile3::C48_mosaic:tile4",
  !  "C48_mosaic:tile3::C48_mosaic:tile5",
  !  "C48_mosaic:tile4::C48_mosaic:tile5",
  !  "C48_mosaic:tile4::C48_mosaic:tile6",
  !  "C48_mosaic:tile5::C48_mosaic:tile6" ;
  ! 
  ! The contact_index contains a pair of four points that defines the two edges that contact to 
  ! each other from the two neighboring tiles.  Assuming the four points are A, B, C, and D.  
  ! A and B defines the edge of tile 1 and C and D defines the edge of tile2.  A is the same point
  ! as C and B is the same as D.  (Ai, Aj) is the index for point A.
  !  Ai:Bi,Aj:Bj::Ci:Di,Cj:Dj
  !
  ! Here is an example of the 12 contacts for a cube-sphere grid
  ! contact_index =
  !  "96:96,1:96::1:1,1:96",
  !  "1:96,96:96::1:1,96:1",
  !  "1:1,1:96::96:1,96:96",
  !  "1:96,1:1::1:96,96:96",
  !  "1:96,96:96::1:96,1:1",
  !  "96:96,1:96::96:1,1:1",
  !  "1:96,1:1::96:96,96:1",
  !  "96:96,1:96::1:1,1:96",
  !  "1:96,96:96::1:1,96:1",
  !  "1:96,96:96::1:96,1:1",
  !  "96:96,1:96::96:1,1:1",
  !  "96:96,1:96::1:1,1:96" ;

#ifdef ESMF_NETCDF
  character(len=128) attstr
  integer   :: i, j, nvars, attlen
  integer   :: tile1, tile2
  integer   :: tiletuple(2,4)
  integer   :: varid1
  integer   :: ncStatus
  character(len=1), allocatable :: contactstring(:,:)
  character(len=1), allocatable :: indexstring(:,:)
  ! Check the attributes first
  ncStatus = nf90_inquire_attribute(ncid, varid, 'standard_name', len=attlen)
  if (CDFCheckError (ncStatus, &
         ESMF_METHOD,  &
         ESMF_SRCLINE, &
         "standard_name attribute does not exit", &
         rc)) return
  ncStatus = nf90_get_att(ncid, varid, 'standard_name', values=attstr)
  if (CDFCheckError (ncStatus, &
         ESMF_METHOD,  &
         ESMF_SRCLINE, &
         "standard_name attribute does not exit", &
         rc)) return
  if (attstr(1:attlen) .ne. 'grid_contact_spec') then
     call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, & 
          msg="- The standard_name attribute is not 'grid_contact_spec'", & 
          ESMF_CONTEXT, rcToReturn=rc) 
     return
  endif
  ! Check the attributes first
  ncStatus = nf90_inquire_attribute(ncid, varid, 'contact_type', len=attlen)
  if (CDFCheckError (ncStatus, &
         ESMF_METHOD,  &
         ESMF_SRCLINE, &
         "contact_type attribute does not exit", &
         rc)) return
  ncStatus = nf90_get_att(ncid, varid, 'contact_type', values=attstr)
  if (CDFCheckError (ncStatus, &
         ESMF_METHOD,  &
         ESMF_SRCLINE, &
         "error reading contact_type attribute", &
         rc)) return
  if (attstr(1:attlen) .ne. 'boundary') then
     call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, & 
          msg="- Only 'boundary' contact_type is supported", & 
          ESMF_CONTEXT, rcToReturn=rc) 
     return
  endif
  allocate(contactstring(dims(1),dims(2)))
  ncStatus = nf90_get_var(ncid, varid, contactstring, start=(/1,1/), count=(/dims(1), dims(2)/))
  if (CDFCheckError (ncStatus, &
         ESMF_METHOD,  &
         ESMF_SRCLINE, &
         "Error reading contacts variable", &
         rc)) return
  do j=1, dims(2)
    call parse_contact(ESMF_UtilArray2String(contactstring(:,j)), tilenames, tile1, tile2)
    contact(1,j)=tile1
    contact(2,j)=tile2
  enddo

  ! read the contact_index
  ncStatus = nf90_inquire_attribute(ncid, varid, 'contact_index', len=attlen)
  if (CDFCheckError (ncStatus, &
         ESMF_METHOD,  &
         ESMF_SRCLINE, &
         "contact_index attribute does not exit", &
         rc)) return
  ncStatus = nf90_get_att(ncid, varid, 'contact_index', values=attstr)
  if (CDFCheckError (ncStatus, &
         ESMF_METHOD,  &
         ESMF_SRCLINE, &
         "error reading contact_index attribute", &
         rc)) return

  ncStatus = nf90_inq_varid(ncid, attstr(1:attlen), varid1)
  if (CDFCheckError (ncStatus, &
         ESMF_METHOD,  &
         ESMF_SRCLINE, &
         attstr(1:attlen), &
         rc)) return
  allocate(indexstring(dims(1),dims(2)))
  ncStatus = nf90_get_var(ncid, varid1, indexstring, start=(/1,1/), count=(/dims(1), dims(2)/))
  if (CDFCheckError (ncStatus, &
         ESMF_METHOD,  &
         ESMF_SRCLINE, &
         "contact_type attribute does not exit", &
         rc)) return
  do j=1, dims(2)
    call parse_contactindex(ESMF_UtilArray2String(indexstring(:,j)), tiletuple)
    connindex(:,:,j)=tiletuple
  enddo

#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 readContacts