PRO  read_all_grid_hdf, gi, field_names, all, VERBOSE = verbose, $
                        STRIDE=stride, START=start, COUNT=count, $
                        particle_density=particle_density
; if file_name is not an array of file names it assume 
; that file_name is a wilcard and read hdf files that match wildcard
; file_name
; with all data matching field_names of each grid 
; and return it in data_structure all
;
  readall        = {readall,                         $
                    file_name: "",                   $
                    Np_total: 0L,                    $
                    dims: [0,0,0],                   $
                    li: [0,0,0],                     $  ; left index
                    ri: [0,0,0],                     $  ; right index
                    N_dfields: 0B,                   $
                    DATA_FIELDS: STRARR(100),        $
                    Data:     PTRARR(100)   }

  initialize, verbose, 0
  files = gi.baryon_file
  nf = N_elements(files)
  if (nf eq 1) then begin
     files = File_Search(files)  
     nf = N_elements(files)
  endif
  all = replicate(readall, nf)

  if (n_elements(count) gt 0) then begin
     if not(max((size(count))[1:2] eq nf)) then begin
        print, "read_all_grid_hdf: count does not have the right dims.", size(count),nf
;        return
    endif
 endif

  if (n_elements(start) gt 0) then begin
     if not(max((size(start))[1:2] eq nf)) then begin
        print, "read_all_grid_hdf: start does not have the right dims."
        return
    endif
endif 

  if (n_elements(stride) gt 0) then begin
     if not(max((size(stride))[1:2] eq nf)) then begin
        print, "read_all_grid_hdf: stride does not have the right dims."
        return
    endif
 endif

  stride_defined = n_elements(stride) gt 0
  start_defined = n_elements(start) gt 0
  count_defined = n_elements(count) gt 0

  if keyword_set(verbose) then print, 'read_all_grid_hdf: reading', field_names, ' from', nf,' files.'
  if files[0] eq '' then BEGIN
      print, 'read_all_grid_hdf: the file "', file_name, '" does not exist!'
      print, ' returning nothing useful'
      RETURN
  ENDIF

  N_df_requested = N_ELEMENTS(field_names)

  special = ['x', 'y', 'z','dx','dy','dz','volume','mass']
  for fi=0L,nf-1 DO BEGIN
     file_name = files[fi]
     sd_id = HDF_SD_START(file_name, /READ)
     cg = gi[fi]
     dd = (cg.right_edge-cg.left_edge)/(cg.end_index-cg.start_index+1)
     ra = cg.end_index-cg.start_index + 1
     d = ra
     this_stride = (stride_defined) ? reform(stride[fi,*]) : $
                   (intarr(n_elements(ra)) + 1)
     this_start = ((start_defined) ? reform(start[fi,*]) : intarr(n_elements(ra)) ) > 0
              

     FOR i=0,N_elements(field_names)-1 DO BEGIN
        field = field_names[i]
        if field eq 'mass' then field = 'Density'
        if field eq 'particle massdensity' then field = 'particle mass'
        if where(field eq special) ne -1 then begin ; this is a special field
           if count_defined then ra = (count[fi,*])
           in = INDGEN(product(ra), /l64)
           data = fltarr(ra[0],ra[1],ra[2])
           ind = array_indices(ra,in, /dimensions)
           case field of
              'x':  data[*] = cg.Left_edge[0] + (0.5D + DOUBLE(this_start[0]+ind[0,*]))*dd[0]
              'y':  data[*] = cg.Left_edge[1] + (0.5D + DOUBLE(this_start[1]+ind[1,*]))*dd[1]
              'z':  data[*] = cg.Left_edge[2] + (0.5D + DOUBLE(this_start[2]+ind[2,*]))*dd[2]
              'dx': data[*] = dd[0]
              'dy': data[*] = dd[1]
              'dz': data[*] = dd[2]
              'volume': data[*] = product(dd)
              ELSE: data[*] = 1.
           endcase
        ENDIF ELSE BEGIN ; regular field to read from disk
           sd_ind = HDF_SD_NAMETOINDEX(sd_id, field)  
           if sd_ind eq -1 then begin
              print, 'read_all_grid_hdf: no such field:', field_names[i], ' in ', file_name
              break
           end else begin
              sds_id = HDF_SD_SELECT(sd_id, sd_ind)
              HDF_SD_GETINFO, sds_id, label = l, dims=d, type=t
              IF ((verbose GT 0) and (fi eq 0)) THEN $
                 if verbose then print, ' read ',l, ' data from file(s) including', file_name
              if N_elements(d) eq 2 then $
                 HDF_SD_GETDATA, sds_id, data, start=this_start[0:1], stride=this_stride[0:1] $
              else begin
                 if (total(this_start lt d) eq 3) then begin
                    if count_defined then begin
                       HDF_SD_GETDATA, sds_id, data, start=this_start, stride=this_stride, $
                                       count=reform(count[fi,*])
                    endif else begin
                       HDF_SD_GETDATA, sds_id, data, start=this_start, stride=this_stride
                    endelse
                 endif else continue
              endelse
              HDF_SD_ENDACCESS, sds_id
           endelse
        ENDELSE
        all[fi].file_name   = file_name
        all[fi].np_total    = N_elements(data)
        s = size(data)
        if (s[0] eq 2) then s[3] = 1
        if (s[0] eq 1) then s[2:3] = 1
        ddd = fltarr(3)  &  ddd[0] = d[*]
        all[fi].li[*]       = this_start[*]
        all[fi].ri[*]       = (s[1:3]+this_start-1)[*] 
        all[fi].dims        = all[fi].ri - all[fi].li + 1
        all[fi].N_Dfields   = N_df_requested
        all[fi].data_fields[i]=field_names[i]
;        print, field_names[i]
        if field_names[i] eq 'particle massdensity' then data *= product(dd)
        if field_names[i] eq 'mass' then data *= product(dd)
        all[fi].Data[i]     = ptr_new(data,/no_copy)
     ENDFOR                     ; loop over field_names
     HDF_SD_END, sd_id
     if (fi mod 10) eq 0 then $
        if check_cancel() then break
     
  ENDFOR                        ; loop over files
;breaki
RETURN
END
  
;  fname = '../../RedshiftOutput0015.dir/Red*grid000*'
;  read_all_grid_hdf, fname, [1,6], data
;END