function find_max_value_gadget
@common_blocks.inc
  densind = where(*parti.names eq 'Density')
  region = where((*parti.d[2]) gt 16 and (*parti.d[2]) lt 18 and $ 
                 (*parti.d[1]) gt 14 and (*parti.d[1]) lt 18 and $ 
                 (*parti.d[0]) gt 14 and (*parti.d[0]) lt 18  $ 
                )
  if densind[0] ne -1 then begin
     dummy = max((*parti.d[densind[0]])[region], ind)
     ind = region[ind]
     xyz_vector = [(*parti.d[0])[ind[0]],(*parti.d[1])[ind[0]],(*parti.d[2])[ind[0]]]
     Print, 'Maximum value is:', dummy, ' at', strcompress(xyz_vector)
  endif
  return, xyz_vector
end
     
pro find_max_value, option,  u_g, field_name, xyz_vector
; using all grid in the array of grid structure: grid_info
; find the spatial coord[0..1,0..1,0..1]  of the highest value of the
; quantitiy specified by sds_num (see hdf files)
; and return it in xyz_vector
  dim_grid_info = size(u_g)
  num_of_grids  = dim_grid_info(1)
  IF (num_of_grids lt 1) THEN BEGIN
      print, 'find_max_value: no grid information in grid_info:', u_g
      STOP
  ENDIF

  data_format = determine_data_format_from_file_name(u_g[0].baryon_file)
  if (data_format eq 2) then begin
     xyz_vector = find_max_value_gadget()
     return
  endif

  print, 'STRING', option

  IF (option eq 'Find max on finest level') THEN BEGIN
      sec_finest_level = max(u_g.level)   - 1 
      grid_info    = u_g(where(u_g.level gt sec_finest_level))
  END ELSE grid_info = u_g

  dim_grid_info = size(grid_info)
  num_of_grids  = dim_grid_info(1)


  a = get_data(grid_info, field_name)
  max_value = -1.e30 ; that's small 
  For i=0L,num_of_grids-1 DO BEGIN
      data_points = grid_info(i).End_index-grid_info(i).Start_index+1
      delta_dist  = (DOUBLE(grid_info(i).Right_edge)- $
                     DOUBLE(grid_info(i).Left_edge))/DOUBLE(data_points)
      max_value_tg =  max(*a[i].data[0],full_index)
      if (max_value_tg ge max_value) THEN BEGIN
          max_value  =  max_value_tg
          index = indgen(3)
          index(0)   =  full_index mod data_points(0)
          index(1)   =  full_index/data_points(0) mod data_points(1)
          index(2)   =  full_index/ $
            (data_points(0) * data_points(1)) mod data_points(2)
          xyz_vector = DOUBLE(grid_info[i].Left_edge) + (index)*delta_dist
          print, 'xyz_vec, max, maxcheck:', xyz_vector, max_value,$
            (*a[i].data[0])[index(0),index(1),index(2)]
      ENDIF
  ENDFOR

  Print, 'Maximum value is:', max_value
END
;.compile find_max_value.pro