def lithospheric_stress(step, trench, ridge, time):
    """calculate stress in the lithosphere"""
    timestep = step.isnap
    base_lith = step.geom.rcmb + 1 - 0.105

    stressfld = step.fields['sII'][0, :, :, 0]
    stressfld = np.ma.masked_where(step.geom.r_mesh[0] < base_lith, stressfld)

    # stress integration in the lithosphere
    dzm = (step.geom.r_coord[1:] - step.geom.r_coord[:-1])
    stress_lith = np.sum((stressfld[:, 1:] * dzm.T), axis=1)
    ph_coord = step.geom.p_coord  # probably doesn't need alias

    # plot stress in the lithosphere
    fig, axis, _, _ = field.plot_scalar(step, 'sII', stressfld,
                                        cmap='plasma_r', 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)
    misc.saveplot(fig, 'lith', timestep)

    # velocity
    vphi = step.fields['v2'][0, :, :, 0]
    vph2 = 0.5 * (vphi + np.roll(vphi, 1, 0))  # interpolate to the same phi

    # position of continents
    concfld = step.fields['c'][0, :, :, 0]
    if step.sdat.par['boundaries']['air_layer']:
        # we are a bit below the surface; delete "-some number"
        # to be just below
        dsa = step.sdat.par['boundaries']['air_thickness']
        # depth to detect the continents
        indcont = np.argmin(abs((1 - dsa) - step.geom.r_coord)) - 10
    else:
        # depth to detect continents
        indcont = -1
    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

    # plot integrated stress in the lithosphere
    fig0, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
    ax1.plot(ph_coord[:-1], vph2[:-1, -1], 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)

    intstr_scale = step.sdat.scales.stress * step.sdat.scales.length / 1.e12
    ax2.plot(ph_coord, stress_lith * intstr_scale, color='k', label='Stress')
    ax2.set_ylabel(r"Integrated stress [$TN\,m^{-1}$]")

    plot_plate_limits(ax1, ridge, trench, conf.plates.vmin,
                      conf.plates.vmax)
    plot_plate_limits(ax2, ridge, trench, conf.plates.stressmin,
                      conf.plates.lstressmax)
    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.stressmin,
        conf.plates.lstressmax, facecolor='#8b6914', alpha=0.2)
    ax2.set_ylim(conf.plates.stressmin, conf.plates.lstressmax)

    misc.saveplot(fig0, 'svelslith', timestep)
calculate stress in the lithosphere