docs: def main_plates(sdat):
      """Plot several plates information."""
      # calculating averaged horizontal surface velocity
      # needed for redimensionalisation
      ilast = sdat.rprof.index.levels[0][-1]
      rlast = sdat.rprof.loc[ilast]
      nprof = 0
      uprof_averaged = rlast.loc[:, 'vhrms'] * 0
      for step in sdat.walk.filter(rprof=True):
          uprof_averaged += step.rprof['vhrms']
          nprof += 1
      uprof_averaged /= nprof
      radius = rlast['r'].values
      if sdat.par['boundaries']['air_layer']:
          dsa = sdat.par['boundaries']['air_thickness']
          isurf = np.argmin(abs(radius - radius[-1] + dsa))
          vrms_surface = uprof_averaged.iloc[isurf]
          isurf = np.argmin(abs((1 - dsa) - radius))
          isurf -= 4  # why different isurf for the rest?
      else:
          isurf = -1
          vrms_surface = uprof_averaged.iloc[isurf]
  
      with misc.InchoateFiles(8, 'plates') as fids:
          fids.fnames = ['plate_velocity', 'distance_subd', 'continents',
                         'flux', 'topography', 'age', 'velderiv', 'velocity']
          fids[0].write('#  it  time  ph_trench vel_trench age_trench\n')
          fids[1].write('#  it      time   time [My]   distance     '
                        'ph_trench     ph_cont  age_trench [My]\n')
  
          istart, iend = None, None
          for step in sdat.walk.filter(fields=['T']):
              # could check other fields too
              timestep = step.isnap
              istart = timestep if istart is None else istart
              iend = timestep
              print('Treating snapshot', timestep)
  
              rcmb = step.geom.rcmb
              # topography
              fname = sdat.filename('sc', timestep=timestep, suffix='.dat')
              topo = np.genfromtxt(str(fname))
              # rescaling topography!
              if sdat.par['boundaries']['air_layer']:
                  topo[:, 1] = topo[:, 1] / (1. - dsa)
  
              time = step.geom.ti_ad * vrms_surface *\
                  conf.scaling.ttransit / conf.scaling.yearins / 1.e6
              trenches, ridges, agetrenches, _, _ =\
                  detect_plates(step, vrms_surface, fids, time)
              plot_plates(step, time, vrms_surface, trenches, ridges,
                          agetrenches, topo, fids)
  
              # prepare for continent plotting
              concfld = step.fields['c'][0, :, :, 0]
              continentsfld = np.ma.masked_where(
                  concfld < 3, concfld)  # plotting continents, to-do
              continentsfld = continentsfld / continentsfld
  
              temp = step.fields['T'][0, :, :, 0]
              tgrad = (temp[:, isurf - 1] - temp[:, isurf]) /\
                  (step.geom.r_coord[isurf] - step.geom.r_coord[isurf - 1])
  
              io_surface(timestep, time, fids[2], concfld[:-1, isurf])
              io_surface(timestep, time, fids[3], tgrad)
              io_surface(timestep, time, fids[4], topo[:, 1])
              if 'age' in conf.plates.plot:
                  io_surface(timestep, time, fids[5],
                             step.fields['age'][0, :, isurf, 0])
  
              # plot viscosity field with position of trenches and ridges
              etamin, _ = sdat.scale(1e-2, 'Pa')
              etamax, _ = sdat.scale(sdat.par['viscosity']['eta_max'], 'Pa')
              fig, axis, _, _ = field.plot_scalar(step, 'eta',
                                                  vmin=etamin, vmax=etamax)
  
              # plotting continents
              field.plot_scalar(step, 'c', continentsfld, axis, False,
                                cmap='cool_r', vmin=0, vmax=0)
              cmap2 = plt.cm.ocean
              cmap2.set_over('m')
  
              # plotting velocity vectors
              field.plot_vec(axis, step, 'v')
  
              # Annotation with time and step
              axis.text(1., 0.9, str(round(time, 0)) + ' My',
                        transform=axis.transAxes)
              axis.text(1., 0.1, str(timestep),
                        transform=axis.transAxes)
  
              # Put arrow where ridges and trenches are
              plot_plate_limits_field(axis, rcmb, ridges, trenches)
  
              misc.saveplot(fig, 'eta', timestep, close=False)
  
              # Zoom
              if conf.plates.zoom is not None:
                  if not 0 <= conf.plates.zoom <= 360:
                      raise error.InvalidZoomError(conf.plates.zoom)
                  if 45 < conf.plates.zoom <= 135:
                      ladd, radd, uadd, dadd = 0.8, 0.8, 0.05, 0.1
                  elif 135 < conf.plates.zoom <= 225:
                      ladd, radd, uadd, dadd = 0.05, 0.1, 0.8, 0.8
                  elif 225 < conf.plates.zoom <= 315:
                      ladd, radd, uadd, dadd = 0.8, 0.8, 0.1, 0.05
                  else:  # >315 or <=45
                      ladd, radd, uadd, dadd = 0.1, 0.05, 0.8, 0.8
                  xzoom = (rcmb + 1) * np.cos(np.radians(conf.plates.zoom))
                  yzoom = (rcmb + 1) * np.sin(np.radians(conf.plates.zoom))
                  axis.set_xlim(xzoom - ladd, xzoom + radd)
                  axis.set_ylim(yzoom - dadd, yzoom + uadd)
                  misc.saveplot(fig, 'etazoom', timestep, close=False)
              plt.close(fig)
  
              # plot stress field with position of trenches and ridges
              if 'str' in conf.plates.plot:
                  fig, axis, _, _ = field.plot_scalar(step, 'sII',
                                                      vmin=0, vmax=300)
  
                  # Annotation with time and step
                  axis.text(1., 0.9, str(round(time, 0)) + ' My',
                            transform=axis.transAxes)
                  axis.text(1., 0.1, str(timestep),
                            transform=axis.transAxes)
  
                  # Put arrow where ridges and trenches are
                  plot_plate_limits_field(axis, rcmb, ridges, trenches)
  
                  misc.saveplot(fig, 's', timestep, close=False)
  
                  # Zoom
                  if conf.plates.zoom is not None:
                      axis.set_xlim(xzoom - ladd, xzoom + radd)
                      axis.set_ylim(yzoom - dadd, yzoom + uadd)
                      misc.saveplot(fig, 'szoom', timestep, close=False)
                  plt.close(fig)
  
                  # calculate stresses in the lithosphere
                  lithospheric_stress(step, trenches, ridges, time)
  
              # plotting the principal deviatoric stress field
              if 'sx' in conf.plates.plot:
                  fig, axis, _, _ = field.plot_scalar(step, 'sII',
                                                      alpha=0.1)
  
                  # plotting continents
                  field.plot_scalar(step, 'c', continentsfld, axis, False,
                                    cmap='cool_r', vmin=0, vmax=0)
                  cmap2 = plt.cm.ocean
                  cmap2.set_over('m')
  
                  # plotting principal deviatoric stress
                  field.plot_vec(axis, step, 'sx')
  
                  # Annotation with time and step
                  axis.text(1., 0.9, str(round(time, 0)) + ' My',
                            transform=axis.transAxes)
                  axis.text(1., 0.1, str(timestep),
                            transform=axis.transAxes)
  
                  # Put arrow where ridges and trenches are
                  plot_plate_limits_field(axis, rcmb, ridges, trenches)
  
                  misc.saveplot(fig, 'sx', timestep)
  
          # determine names of files
          ptn = misc.out_name('{}_{}_{}')
          stem = ptn.format(fids.fnames[0], istart, iend)
          idx = 0
          fmt = '{}.dat'
          while pathlib.Path(fmt.format(stem, idx)).is_file():
              fmt = '{}_{}.dat'
              idx += 1
          fids.fnames = [fmt.format(ptn.format(fname, istart, iend), idx)
                         for fname in fids.fnames]