docs: def t_surfabu(self, num_frame, xax, t0_model=0,
                    title='surface abundance', t_eps=1.e-3,
                    plot_CO_ratio=False):
          """
          t_surfabu plots surface abundance evolution as a function of
          time.
  
          Parameters
          ----------
          num_frame : integer
              Number of frame to plot this plot into, if <0 don't open
              figure.
          xax : string
              Either model, time or logrevtime to indicate what is to be
              used on the x-axis.
          t0_model : integer, optional
              Model for the zero point in time, for AGB plots this would
              be usually the model of the 1st TP, which can be found with
              the Kippenhahn plot.  The default is 0.
          title : string, optional
              Figure title.  The default is "surface abundance".
          t_eps : float, optional
              Time eps at end for logrevtime.  The default is 1.e-3.
          plot_CO_ratio : boolean, optional
              On second axis True/False.  The default is False.
  
          """
          if num_frame >= 0:
              pyl.figure(num_frame)
  
          if xax == 'time':
              xaxisarray = self.get('star_age')[t0_model:]
          elif xax == 'model':
              xaxisarray = self.get('model_number')[t0_model:]
          elif xax == 'logrevtime':
              xaxisarray = self.get('star_age')
              xaxisarray=np.log10(max(xaxisarray[t0_model:])+t_eps-xaxisarray[t0_model:])
          else:
              print('t-surfabu error: invalid string for x-axis selction.'+ \
                    ' needs to be "time" or "model"')
  
          star_mass         = self.get('star_mass')
          surface_c12       = self.get('surface_c12')
          surface_c13       = self.get('surface_c13')
          surface_n14       = self.get('surface_n14')
          surface_o16       = self.get('surface_o16')
  
          target_n14 = -3.5
  
  
          COratio=old_div((surface_c12*4.),(surface_o16*3.))
          t0_mod=xaxisarray[t0_model]
          log10_c12=np.log10(surface_c12[t0_model:])
  
          symbs=['k:','-','--','-.','b:','-','--','k-.',':','-','--','-.']
  
          pyl.plot(xaxisarray,log10_c12,\
                       symbs[0],label='$^{12}\mathrm{C}$')
          pyl.plot(xaxisarray,np.log10(surface_c13[t0_model:]),\
                       symbs[1],label='$^{13}\mathrm{C}$')
          pyl.plot(xaxisarray,np.log10(surface_n14[t0_model:]),\
                       symbs[2],label='$^{14}\mathrm{N}$')
          pyl.plot(xaxisarray,np.log10(surface_o16[t0_model:]),\
                       symbs[3],label='$^{16}\mathrm{O}$')
  #                pyl.plot([min(xaxisarray[t0_model:]-t0_mod),max(xaxisarray[t0_model:]-t0_mod)],[target_n14,target_n14])
  
          pyl.ylabel('mass fraction $\log X$')
          pyl.legend(loc=2)
  
          if xax == 'time':
              pyl.xlabel('t / yrs')
          elif xax == 'model':
              pyl.xlabel('model number')
          elif xax == 'logrevtime':
              pyl.xlabel('$\\log t-tfinal$')
          if plot_CO_ratio:
              pyl.twinx()
              pyl.plot(xaxisarray,COratio[t0_model:],'-k',label='CO ratio')
              pyl.ylabel('C/O ratio')
              pyl.legend(loc=4)
          pyl.title(title)
          if xax == 'logrevtime':
              self._xlimrev()