docs: def plot_scalar(step, var, field=None, axis=None, set_cbar=True, **extra):
      """Plot scalar field.
  
      Args:
          step (:class:`~stagpy.stagyydata._Step`): a step of a StagyyData
              instance.
          var (str): the scalar field name.
          field (:class:`numpy.array`): if not None, it is plotted instead of
              step.fields[var].  This is useful to plot a masked or rescaled
              array.  Note that if conf.scaling.dimensional is True, this
              field will be scaled accordingly.
          axis (:class:`matplotlib.axes.Axes`): the axis objet where the field
              should be plotted.  If set to None, a new figure with one subplot
              is created.
          set_cbar (bool): whether to add a colorbar to the plot.
          extra (dict): options that will be passed on to
              :func:`matplotlib.axes.Axes.pcolormesh`.
      Returns:
          fig, axis, surf, cbar
              handles to various :mod:`matplotlib` objects, respectively the
              figure, the axis, the surface returned by
              :func:`~matplotlib.axes.Axes.pcolormesh`, and the colorbar returned
              by :func:`matplotlib.pyplot.colorbar`.
      """
      if var in phyvars.FIELD:
          meta = phyvars.FIELD[var]
      else:
          meta = phyvars.FIELD_EXTRA[var]
          meta = phyvars.Varf(misc.baredoc(meta.description), meta.dim)
      if step.geom.threed:
          raise NotAvailableError('plot_scalar only implemented for 2D fields')
  
      xmesh, ymesh, fld = get_meshes_fld(step, var)
      xmin, xmax = xmesh.min(), xmesh.max()
      ymin, ymax = ymesh.min(), ymesh.max()
  
      if field is not None:
          fld = field
      if conf.field.perturbation:
          fld = fld - np.mean(fld, axis=0)
      if conf.field.shift:
          fld = np.roll(fld, conf.field.shift, axis=0)
  
      fld, unit = step.sdat.scale(fld, meta.dim)
  
      if axis is None:
          fig, axis = plt.subplots(ncols=1)
      else:
          fig = axis.get_figure()
  
      if step.sdat.par['magma_oceans_in']['magma_oceans_mode']:
          rcmb = step.sdat.par['geometry']['r_cmb']
          xmax = rcmb + 1
          ymax = xmax
          xmin = -xmax
          ymin = -ymax
          rsurf = xmax if step.timeinfo['thick_tmo'] > 0 \
              else step.geom.r_mesh[0, 0, -3]
          cmb = mpat.Circle((0, 0), rcmb, color='dimgray', zorder=0)
          psurf = mpat.Circle((0, 0), rsurf, color='indianred', zorder=0)
          axis.add_patch(psurf)
          axis.add_patch(cmb)
  
      extra_opts = dict(
          cmap=conf.field.cmap.get(var),
          vmin=conf.plot.vmin,
          vmax=conf.plot.vmax,
          norm=mpl.colors.LogNorm() if var == 'eta' else None,
          rasterized=conf.plot.raster,
          shading='gouraud' if conf.field.interpolate else 'flat',
      )
      extra_opts.update(extra)
      surf = axis.pcolormesh(xmesh, ymesh, fld, **extra_opts)
  
      cbar = None
      if set_cbar:
          cbar = plt.colorbar(surf, shrink=conf.field.shrinkcb)
          cbar.set_label(meta.description +
                         (' pert.' if conf.field.perturbation else '') +
                         (' ({})'.format(unit) if unit else ''))
      if step.geom.spherical or conf.plot.ratio is None:
          plt.axis('equal')
          plt.axis('off')
      else:
          axis.set_aspect(conf.plot.ratio / axis.get_data_ratio())
      axis.set_adjustable('box')
      axis.set_xlim(xmin, xmax)
      axis.set_ylim(ymin, ymax)
      return fig, axis, surf, cbar