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],
        ))