def lithospheric_stress(step, trench, ridge, time):
timestep = step.isnap<EOL>base_lith = step.geom.rcmb + <NUM_LIT:1> - <NUM_LIT><EOL>stressfld = step.fields['<STR_LIT>'][<NUM_LIT:0>, :, :, <NUM_LIT:0>]<EOL>stressfld = np.ma.masked_where(step.geom.r_mesh[<NUM_LIT:0>] < base_lith, stressfld)<EOL>dzm = (step.geom.r_coord[<NUM_LIT:1>:] - step.geom.r_coord[:-<NUM_LIT:1>])<EOL>stress_lith = np.sum((stressfld[:, <NUM_LIT:1>:] * dzm.T), axis=<NUM_LIT:1>)<EOL>ph_coord = step.geom.p_coord  <EOL>fig, axis, _, _ = field.plot_scalar(step, '<STR_LIT>', stressfld,<EOL>cmap='<STR_LIT>', vmin=<NUM_LIT:0>, vmax=<NUM_LIT>)<EOL>axis.text(<NUM_LIT:1.>, <NUM_LIT>, str(round(time, <NUM_LIT:0>)) + '<STR_LIT>', transform=axis.transAxes)<EOL>axis.text(<NUM_LIT:1.>, <NUM_LIT:0.1>, str(timestep), transform=axis.transAxes)<EOL>misc.saveplot(fig, '<STR_LIT>', timestep)<EOL>vphi = step.fields['<STR_LIT>'][<NUM_LIT:0>, :, :, <NUM_LIT:0>]<EOL>vph2 = <NUM_LIT:0.5> * (vphi + np.roll(vphi, <NUM_LIT:1>, <NUM_LIT:0>))  <EOL>concfld = step.fields['<STR_LIT:c>'][<NUM_LIT:0>, :, :, <NUM_LIT:0>]<EOL>if step.sdat.par['<STR_LIT>']['<STR_LIT>']:<EOL><INDENT>dsa = step.sdat.par['<STR_LIT>']['<STR_LIT>']<EOL>indcont = np.argmin(abs((<NUM_LIT:1> - dsa) - step.geom.r_coord)) - <NUM_LIT:10><EOL><DEDENT>else:<EOL><INDENT>indcont = -<NUM_LIT:1><EOL><DEDENT>if step.sdat.par['<STR_LIT>']['<STR_LIT>'] andnot step.sdat.par['<STR_LIT>']['<STR_LIT>']:<EOL><INDENT>continents = np.ma.masked_where(<EOL>np.logical_or(concfld[:-<NUM_LIT:1>, indcont] < <NUM_LIT:3>,<EOL>concfld[:-<NUM_LIT:1>, indcont] > <NUM_LIT:4>),<EOL>concfld[:-<NUM_LIT:1>, indcont])<EOL><DEDENT>elif step.sdat.par['<STR_LIT>']['<STR_LIT>'] andstep.sdat.par['<STR_LIT>']['<STR_LIT>']:<EOL><INDENT>continents = np.ma.masked_where(<EOL>np.logical_or(concfld[:-<NUM_LIT:1>, indcont] < <NUM_LIT:3>,<EOL>concfld[:-<NUM_LIT:1>, indcont] > <NUM_LIT:5>),<EOL>concfld[:-<NUM_LIT:1>, indcont])<EOL><DEDENT>elif step.sdat.par['<STR_LIT>']['<STR_LIT>']:<EOL><INDENT>continents = np.ma.masked_where(<EOL>concfld[:-<NUM_LIT:1>, indcont] < <NUM_LIT:3>, concfld[:-<NUM_LIT:1>, indcont])<EOL><DEDENT>else:<EOL><INDENT>continents = np.ma.masked_where(<EOL>concfld[:-<NUM_LIT:1>, indcont] < <NUM_LIT:2>, concfld[:-<NUM_LIT:1>, indcont])<EOL><DEDENT>continentsall = continents / continents<EOL>fig0, (ax1, ax2) = plt.subplots(<NUM_LIT:2>, <NUM_LIT:1>, sharex=True, figsize=(<NUM_LIT:12>, <NUM_LIT:8>))<EOL>ax1.plot(ph_coord[:-<NUM_LIT:1>], vph2[:-<NUM_LIT:1>, -<NUM_LIT:1>], label='<STR_LIT>')<EOL>ax1.axhline(y=<NUM_LIT:0>, xmin=<NUM_LIT:0>, xmax=<NUM_LIT:2> * np.pi,<EOL>color='<STR_LIT>', ls='<STR_LIT>', alpha=<NUM_LIT>)<EOL>ax1.set_ylabel("<STR_LIT>")<EOL>ax1.text(<NUM_LIT>, <NUM_LIT>, str(round(time, <NUM_LIT:0>)) + '<STR_LIT>',<EOL>transform=ax1.transAxes)<EOL>ax1.text(<NUM_LIT>, <NUM_LIT>, str(round(step.geom.ti_ad, <NUM_LIT:8>)),<EOL>transform=ax1.transAxes)<EOL>intstr_scale = step.sdat.scales.stress * step.sdat.scales.length / <NUM_LIT><EOL>ax2.plot(ph_coord, stress_lith * intstr_scale, color='<STR_LIT:k>', label='<STR_LIT>')<EOL>ax2.set_ylabel(r"<STR_LIT>")<EOL>plot_plate_limits(ax1, ridge, trench, conf.plates.vmin,<EOL>conf.plates.vmax)<EOL>plot_plate_limits(ax2, ridge, trench, conf.plates.stressmin,<EOL>conf.plates.lstressmax)<EOL>ax1.set_xlim(<NUM_LIT:0>, <NUM_LIT:2> * np.pi)<EOL>ax1.set_title(timestep)<EOL>ax1.fill_between(<EOL>ph_coord[:-<NUM_LIT:1>], continentsall * conf.plates.vmin,<EOL>conf.plates.vmax, facecolor='<STR_LIT>', alpha=<NUM_LIT>)<EOL>ax1.set_ylim(conf.plates.vmin, conf.plates.vmax)<EOL>ax2.fill_between(<EOL>ph_coord[:-<NUM_LIT:1>], continentsall * conf.plates.stressmin,<EOL>conf.plates.lstressmax, facecolor='<STR_LIT>', alpha=<NUM_LIT>)<EOL>ax2.set_ylim(conf.plates.stressmin, conf.plates.lstressmax)<EOL>misc.saveplot(fig0, '<STR_LIT>', timestep)<EOL>
calculate stress in the lithosphere