docs: def read_field_h5(xdmf_file, fieldname, snapshot, header=None):
      """Extract field data from hdf5 files.
  
      Args:
          xdmf_file (:class:`pathlib.Path`): path of the xdmf file.
          fieldname (str): name of field to extract.
          snapshot (int): snapshot number.
          header (dict): geometry information.
      Returns:
          (dict, numpy.array): geometry information and field data. None
              is returned if data is unavailable.
      """
      if header is None:
          header, xdmf_root = read_geom_h5(xdmf_file, snapshot)
      else:
          xdmf_root = xmlET.parse(str(xdmf_file)).getroot()
  
      npc = header['nts'] // header['ncs']  # number of grid point per node
      flds = np.zeros(_flds_shape(fieldname, header))
      data_found = False
  
      for elt_subdomain in xdmf_root[0][0][snapshot].findall('Grid'):
          ibk = int(elt_subdomain.get('Name').startswith('meshYang'))
          for data_attr in elt_subdomain.findall('Attribute'):
              if data_attr.get('Name') != fieldname:
                  continue
              icore, fld = _get_field(xdmf_file, data_attr.find('DataItem'))
              # for some reason, the field is transposed
              fld = fld.T
              shp = fld.shape
              if shp[-1] == 1 and header['nts'][0] == 1:  # YZ
                  fld = fld.reshape((shp[0], 1, shp[1], shp[2]))
                  if header['rcmb'] < 0:
                      fld = fld[(2, 0, 1), ...]
              elif shp[-1] == 1:  # XZ
                  fld = fld.reshape((shp[0], shp[1], 1, shp[2]))
                  if header['rcmb'] < 0:
                      fld = fld[(0, 2, 1), ...]
              elif header['nts'][1] == 1:  # cart XZ
                  fld = fld.reshape((1, shp[0], 1, shp[1]))
              ifs = [icore // np.prod(header['ncs'][:i]) % header['ncs'][i] *
                     npc[i] for i in range(3)]
              if header['zp']:  # remove top row
                  fld = fld[:, :, :, :-1]
              flds[:,
                   ifs[0]:ifs[0] + npc[0] + header['xp'],
                   ifs[1]:ifs[1] + npc[1] + header['yp'],
                   ifs[2]:ifs[2] + npc[2],
                   ibk] = fld
              data_found = True
  
      flds = _post_read_flds(flds, header)
  
      return (header, flds) if data_found else None