     double mass_, posX, posY, velX, velY;
     double posX_old = 0, posY_old;
     double deltaT = p.dt*p.itstp;
-    os << "#Time posX posY velX velY mass Ue Ui Uphi Etot\n";
+    os << "#Time(1) posX(2) posY(3) velX(4) velY(5) mass(6) diff(7) (m_tot-m_0)/m_0(8) "
+       << "Ue(9) Ui(10) Uphi(11) Utot(12) (U_tot-U_0)/U_0(13) diss(14) \n";
     for( unsigned idx=1; idx<=num_out; idx++)
     {
         t5file.get_field( input_h, "electrons", idx);

        double Ue = dg::blas2::dot( input0, w2d, ln0);
        double Ui = p.tau*dg::blas2::dot( input1, w2d, ln1);
        double Uphi = energy[(idx-1)*p.itstp] - Ue - Ui;
        //mass_ = dg::blas2::dot( one, w2d, input0 ); 
        //double mass_p=mass[(idx-1)*p.itstp];
        //if( mass_ - mass_p> 1e-14) 
        //std::cerr<< "Diff masses: "<<mass_ - mass_p<<"\n";
        if( p.global)
            thrust::transform( input0.begin(), input0.end(), input0.begin(), dg::PLUS<double>(-1));
        mass_ = dg::blas2::dot( one, w2d, input0 ); 
        os << t5file.get_time( idx)<<" ";

        posX = dg::blas2::dot( xvec, w2d, input0)/mass_ - p.posX*p.lx;
        posY = dg::blas2::dot( yvec, w2d, input0)/mass_ - p.posY*p.ly;

        velY = (posY - posY_old)/deltaT;
        posX_old = posX;
        posY_old = posY;
        //output
        os << t5file.get_time( idx);//(1)
        os << " "<<posX << " " << posY << " "<<velX<<" "<<velY;//(2-5)
        os << " "<<mass[(idx-1)*p.itstp] << " "<<diffusion[(idx-1)*p.itstp];//(6,7)
        os << " "<< (mass[(idx-1)*p.itstp]-mass[0])/(mass[0]-grid.lx()*grid.ly());//blob mass is mass[] - Area (8)
        os << " "<<Ue<<" "<<Ui<<" "<<Uphi<<" "<<energy[(idx-1)*p.itstp]; //(9-12)
        os << " "<<(energy[(idx-1)*p.itstp]-energy[0])/energy[0];//(13)
        os << " "<<dissipation[(idx-1)*p.itstp]; //(14)
        os <<"\n";
    }
    os.close();
