docs: def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
                  topo, fids):
      """handle ploting stuff"""
      vphi = step.fields['v2'][0, :, :, 0]
      tempfld = step.fields['T'][0, :, :, 0]
      concfld = step.fields['c'][0, :, :, 0]
      timestep = step.isnap
  
      if step.sdat.par['boundaries']['air_layer']:
          dsa = step.sdat.par['boundaries']['air_thickness']
          # we are a bit below the surface; delete "-some number"
          # to be just below
          # the surface (that is considered plane here); should check if you are
          # in the thermal boundary layer
          indsurf = np.argmin(abs((1 - dsa) - step.geom.r_coord)) - 4
          # depth to detect the continents
          indcont = np.argmin(abs((1 - dsa) - np.array(step.geom.r_coord))) - 10
      else:
          indsurf = -1
          indcont = -1  # depth to detect continents
  
      if step.sdat.par['boundaries']['air_layer'] and\
         not step.sdat.par['continents']['proterozoic_belts']:
          continents = np.ma.masked_where(
              np.logical_or(concfld[:-1, indcont] < 3,
                            concfld[:-1, indcont] > 4),
              concfld[:-1, indcont])
      elif (step.sdat.par['boundaries']['air_layer'] and
            step.sdat.par['continents']['proterozoic_belts']):
          continents = np.ma.masked_where(
              np.logical_or(concfld[:-1, indcont] < 3,
                            concfld[:-1, indcont] > 5),
              concfld[:-1, indcont])
      elif step.sdat.par['tracersin']['tracers_weakcrust']:
          continents = np.ma.masked_where(
              concfld[:-1, indcont] < 3, concfld[:-1, indcont])
      else:
          continents = np.ma.masked_where(
              concfld[:-1, indcont] < 2, concfld[:-1, indcont])
  
      # masked array, only continents are true
      continentsall = continents / continents
  
      ph_coord = step.geom.p_coord
  
      # velocity
      vph2 = 0.5 * (vphi + np.roll(vphi, 1, 0))  # interpolate to the same phi
      dvph2 = (np.diff(vph2[:, indsurf]) / (ph_coord[0] * 2.))
  
      # plotting
      fig0, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True, figsize=(12, 8))
      ax1.plot(ph_coord[:-1], concfld[:-1, indsurf], color='g', label='Conc')
      ax2.plot(ph_coord[:-1], tempfld[:-1, indsurf], color='k', label='Temp')
      ax3.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
  
      ax1.fill_between(
          ph_coord[:-1], continents, 1., facecolor='#8B6914', alpha=0.2)
      ax2.fill_between(
          ph_coord[:-1], continentsall, 0., facecolor='#8B6914', alpha=0.2)
  
      tempmin = step.sdat.par['boundaries']['topT_val'] * 0.9\
          if step.sdat.par['boundaries']['topT_mode'] == 'iso' else 0.0
      tempmax = step.sdat.par['boundaries']['botT_val'] * 0.35\
          if step.sdat.par['boundaries']['botT_mode'] == 'iso' else 0.8
  
      ax2.set_ylim(tempmin, tempmax)
      ax3.fill_between(
          ph_coord[:-1], continentsall * round(1.5 * np.amax(dvph2), 1),
          round(np.amin(dvph2) * 1.1, 1), facecolor='#8B6914', alpha=0.2)
      ax3.set_ylim(conf.plates.vmin, conf.plates.vmax)
  
      ax1.set_ylabel("Concentration")
      ax2.set_ylabel("Temperature")
      ax3.set_ylabel("Velocity")
      ax1.set_title(timestep)
      ax1.text(0.95, 1.07, str(round(time, 0)) + ' My',
               transform=ax1.transAxes)
      ax1.text(0.01, 1.07, str(round(step.geom.ti_ad, 8)),
               transform=ax1.transAxes)
  
      plot_plate_limits(ax3, ridge, trench, conf.plates.vmin,
                        conf.plates.vmax)
  
      misc.saveplot(fig0, 'sveltempconc', timestep)
  
      # plotting velocity and velocity derivative
      fig0, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
      ax1.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
      ax1.axhline(y=0, xmin=0, xmax=2 * np.pi,
                  color='black', ls='solid', alpha=0.2)
      ax1.set_ylabel("Velocity")
      ax1.text(0.95, 1.07, str(round(time, 0)) + ' My',
               transform=ax1.transAxes)
      ax1.text(0.01, 1.07, str(round(step.geom.ti_ad, 8)),
               transform=ax1.transAxes)
      ax2.plot(ph_coord[:-1] + ph_coord[0], dvph2, color='k', label='dv')
      ax2.set_ylabel("dv")
  
      plot_plate_limits(ax1, ridge, trench, conf.plates.vmin,
                        conf.plates.vmax)
      plot_plate_limits(ax2, ridge, trench, conf.plates.dvmin,
                        conf.plates.dvmax)
      ax1.set_xlim(0, 2 * np.pi)
      ax1.set_title(timestep)
  
      ax1.fill_between(
          ph_coord[:-1], continentsall * conf.plates.vmin, conf.plates.vmax,
          facecolor='#8b6914', alpha=0.2)
      ax1.set_ylim(conf.plates.vmin, conf.plates.vmax)
      ax2.fill_between(
          ph_coord[:-1], continentsall * conf.plates.dvmin,
          conf.plates.dvmax, facecolor='#8b6914', alpha=0.2)
      ax2.set_ylim(conf.plates.dvmin, conf.plates.dvmax)
  
      misc.saveplot(fig0, 'sveldvel', timestep)
  
      # plotting velocity and second invariant of stress
      if 'str' in conf.plates.plot:
          stressfld = step.fields['sII'][0, :, :, 0]
          fig0, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
          ax1.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
          ax1.axhline(y=0, xmin=0, xmax=2 * np.pi,
                      color='black', ls='solid', alpha=0.2)
          ax1.set_ylabel("Velocity")
          ax1.text(0.95, 1.07, str(round(time, 0)) + ' My',
                   transform=ax1.transAxes)
          ax1.text(0.01, 1.07, str(round(step.geom.ti_ad, 8)),
                   transform=ax1.transAxes)
          ax2.plot(ph_coord[:-1],
                   stressfld[:-1, indsurf] * step.sdat.scales.stress / 1.e6,
                   color='k', label='Stress')
          ax2.set_ylim(conf.plates.stressmin, conf.plates.stressmax)
          ax2.set_ylabel("Stress [MPa]")
  
          plot_plate_limits(ax1, ridge, trench,
                            conf.plates.vmin, conf.plates.vmax)
          plot_plate_limits(ax2, ridge, trench,
                            conf.plates.stressmin, conf.plates.stressmax)
          ax1.set_xlim(0, 2 * np.pi)
          ax1.set_title(timestep)
  
          ax1.fill_between(
              ph_coord[:-1], continentsall * conf.plates.vmin,
              conf.plates.vmax, facecolor='#8B6914', alpha=0.2)
          ax1.set_ylim(conf.plates.vmin, conf.plates.vmax)
          ax2.fill_between(
              ph_coord[:-1], continentsall * conf.plates.dvmin,
              conf.plates.dvmax,
              facecolor='#8B6914', alpha=0.2)
  
          misc.saveplot(fig0, 'svelstress', timestep)
  
      # plotting velocity
      fig1, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
      ax1.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
      ax1.axhline(y=0, xmin=0, xmax=2 * np.pi,
                  color='black', ls='solid', alpha=0.2)
      ax1.set_ylim(conf.plates.vmin, conf.plates.vmax)
      ax1.set_ylabel("Velocity")
      ax1.text(0.95, 1.07, str(round(time, 0)) + ' My',
               transform=ax1.transAxes)
      plot_plate_limits(ax1, ridge, trench, conf.plates.vmin,
                        conf.plates.vmax)
  
      # plotting velocity and age at surface
      if 'age' in conf.plates.plot:
          agefld = step.fields['age'][0, :, :, 0]
          age_surface = np.ma.masked_where(
              agefld[:, indsurf] < 0.00001, agefld[:, indsurf])
          age_surface_dim = (age_surface * vrms_surface * conf.scaling.ttransit /
                             conf.scaling.yearins / 1.e6)
  
          fig2, (ax3, ax4) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
          ax3.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
          ax3.axhline(
              y=0, xmin=0, xmax=2 * np.pi,
              color='black', ls='solid', alpha=0.2)
          ax3.set_ylim(conf.plates.vmin, conf.plates.vmax)
          ax3.set_ylabel("Velocity")
          ax3.text(0.95, 1.07, str(round(time, 0)) + ' My',
                   transform=ax3.transAxes)
          ax3.fill_between(
              ph_coord[:-1], continentsall * conf.plates.vmax,
              conf.plates.vmin, facecolor='#8B6914', alpha=0.2)
          plot_plate_limits(ax3, ridge, trench,
                            conf.plates.vmin, conf.plates.vmax)
  
      times_subd = []
      age_subd = []
      distance_subd = []
      ph_trench_subd = []
      ph_cont_subd = []
      if step.sdat.par['switches']['cont_tracers']:
          for i, trench_i in enumerate(trench):
              # detection of the distance in between subduction and continent
              ph_coord_noend = ph_coord[:-1]
              angdistance1 = abs(ph_coord_noend[continentsall == 1] - trench_i)
              angdistance2 = 2. * np.pi - angdistance1
              angdistance = np.minimum(angdistance1, angdistance2)
              distancecont = min(angdistance)
              argdistancecont = np.argmin(angdistance)
              continentpos = ph_coord_noend[continentsall == 1][argdistancecont]
  
              ph_trench_subd.append(trench_i)
              age_subd.append(agetrench[i])
              ph_cont_subd.append(continentpos)
              distance_subd.append(distancecont)
              times_subd.append(step.geom.ti_ad)
  
              if angdistance1[argdistancecont] < angdistance2[argdistancecont]:
                  if continentpos - trench_i < 0:  # continent is on the left
                      distancecont = - distancecont
                  ax1.annotate('', xy=(trench_i + distancecont, 2000),
                               xycoords='data', xytext=(trench_i, 2000),
                               textcoords='data',
                               arrowprops=dict(arrowstyle="->", lw="2",
                                               shrinkA=0, shrinkB=0))
              else:  # distance over boundary
                  xy_anot, xy_text = 0, 2 * np.pi
                  if continentpos - trench_i < 0:
                      xy_anot, xy_text = xy_text, xy_anot
                  ax1.annotate('', xy=(xy_anot, 2000),
                               xycoords='data', xytext=(trench_i, 2000),
                               textcoords='data',
                               arrowprops=dict(arrowstyle="-", lw="2",
                                               shrinkA=0, shrinkB=0))
                  ax1.annotate('', xy=(continentpos, 2000),
                               xycoords='data', xytext=(xy_text, 2000),
                               textcoords='data',
                               arrowprops=dict(arrowstyle="->", lw="2",
                                               shrinkA=0, shrinkB=0))
  
      ax1.fill_between(
          ph_coord[:-1], continentsall * conf.plates.vmin,
          conf.plates.vmax, facecolor='#8B6914', alpha=0.2)
      ax2.set_ylabel("Topography [km]")
      ax2.axhline(y=0, xmin=0, xmax=2 * np.pi,
                  color='black', ls='solid', alpha=0.2)
      ax2.plot(topo[:, 0],
               topo[:, 1] * step.sdat.scales.length / 1.e3,
               color='black')
      ax2.set_xlim(0, 2 * np.pi)
      ax2.set_ylim(conf.plates.topomin, conf.plates.topomax)
      ax2.fill_between(
          ph_coord[:-1], continentsall * conf.plates.topomax,
          conf.plates.topomin, facecolor='#8B6914', alpha=0.2)
      plot_plate_limits(ax2, ridge, trench, conf.plates.topomin,
                        conf.plates.topomax)
      ax1.set_title(timestep)
      misc.saveplot(fig1, 'sveltopo', timestep)
  
      if 'age' in conf.plates.plot:
          ax4.set_ylabel("Seafloor age [My]")
          # in dimensions
          ax4.plot(ph_coord[:-1], age_surface_dim[:-1], color='black')
          ax4.set_xlim(0, 2 * np.pi)
          ax4.fill_between(
              ph_coord[:-1], continentsall * conf.plates.agemax,
              conf.plates.agemin, facecolor='#8B6914', alpha=0.2)
          ax4.set_ylim(conf.plates.agemin, conf.plates.agemax)
          plot_plate_limits(ax4, ridge, trench, conf.plates.agemin,
                            conf.plates.agemax)
          ax3.set_title(timestep)
          misc.saveplot(fig2, 'svelage', timestep)
  
      # writing the output into a file, all time steps are in one file
      for isubd in np.arange(len(distance_subd)):
          fids[1].write("%6.0f %11.7f %11.3f %10.6f %10.6f %10.6f %11.3f\n" % (
              timestep,
              times_subd[isubd],
              time,
              distance_subd[isubd],
              ph_trench_subd[isubd],
              ph_cont_subd[isubd],
              age_subd[isubd],
          ))