<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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()