docs: 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)