function grid_int_coords, INFO=info
    ; start with the root
    ; 64 levels of refinement should be enough for anybody.
    dims = LON64ARR(6,N_ELEMENTS(info))
    l_max = max(info.level)
    h = histogram([info.level], REVERSE_INDICES=ri)
    num = info.num
    for l=0,l_max do begin
        level_indices = ri[ri[l]:ri[l+1]-1]
        for ii=0,n_elements(level_indices)-1 do begin
            i = level_indices[ii]
            g = info[i]
            p = (where(num eq g.parent))[0]
            if p lt 0 then $
                dims[0,i] = [0,0,0,g.dim-1] $
            else begin
                scale = ishft(1L,g.level-info[p].level)
                left = ISHFT(g.parent_s_index + dims[0:2, p], g.level-info[p].level)
                dims[0,i] = [left, left + g.dim - 1]
            endelse
        endfor
    endfor
    return, dims
end

function grid_neighbors, INFO=info, DIMS=dims, CELLS=cells
    l_max = max(info.level)
    h = histogram([info.level], REVERSE_INDICES=ri)
    max_neighbors = max(h)
    neighbors = LONARR(max_neighbors, N_ELEMENTS(info))
    neighbors[*] = -1
    ; for each grid, find the list of possible neighbors
    for k=0,l_max do begin
        level_indices = ri[ri[k]:ri[k+1]-1]
        level_info = info[level_indices]
        level_dims = dims[*,level_indices]
        for i=0L,N_ELEMENTS(level_indices)-1 do begin
            l = CMREPLICATE(level_dims[0:2,i], N_ELEMENTS(level_indices)) - cells
            r = CMREPLICATE(level_dims[3:5,i], N_ELEMENTS(level_indices)) + cells
            ; the expression (b-c)*(d-a) will be positive iff ranges [a,b] and [c,d] overlap
            ; we need three overlaps and/or contacts for a grid to be adjacent
            ind = where((TOTAL((r-level_dims[0:2,*])*(level_dims[3:5,*]-l) ge 0L, 1, /INTEGER) eq 3) and (level_info.num ne level_info[i].num))
            if ind[0] ge 0 then $
                neighbors[0, level_indices[i]] = level_indices[ind]
        endfor
    endfor
    return, neighbors
end