     dg::blas1::axpby( -1., omega, 1., chi); //chi  = + N_i U_i -N_e U_e
     dg::blas1::axpby( -1., y[2], 1., y[3], omega); //omega  = - U_e + U_i   
     double Dres = -p.c*dg::blas2::dot(omega, w3d, chi); //- C*(N_i U_i + N_e U_e)(U_i - U_e)
-    double nu_parallel[] = {-p.mu[0]/p.c, -p.mu[0]*p.tau[1]/p.c, p.nu_parallel, p.nu_parallel};
-    for( unsigned i=0; i<2;i++)
-    {
-        //Compute parallel dissipative energy for N/////////////////////////////
-        if( p.nu_parallel != 0)
-        {
-            if (p.pollim==1) dsN_.set_boundaries( p.bc, 0, 0);             
-            if (p.pardiss==0)
-            {
-                dsN_.forward(y[i], omega); 
-                dsN_.forwardT(omega,lambda);
-                dg::blas1::axpby( 0.5*nu_parallel[i], lambda, 0.,lambda,lambda); 
-
-                dsN_.backward( y[i], omega); 
-                dsN_.backwardT(omega,lambda);
-                dg::blas1::axpby( 0.5*nu_parallel[i], lambda,  1., lambda,lambda); 
-            }
-            if (p.pardiss==1)
-            {
-//                     dsN_.dsz(y[i],omega);                                            //ds^2 N 
-//                     dg::blas1::axpby( p.nu_parallel, omega, 0., lambda,lambda);     //lambda = nu_para*ds^2 N 
-//                     dsN_(y[i], dsy[i]);       
-//                     dg::blas1::pointwiseDot(gradlnB, dsy[i], omega);                // ds lnB ds N
-//                     dg::blas1::axpby(-p.nu_parallel, omega, 1., lambda,lambda);     // lambda += nu_para*ds lnB ds N
-                dsN_.forward( y[i], omega); 
-                dsDIR_.forwardTD(omega,lambda);
-                dg::blas1::axpby( 0.5*nu_parallel[i], lambda, 1., lambda,lambda); 
-                dsN_.backward( y[i], omega); 
-                dsDIR_.backwardTD(omega,lambda);
-                dg::blas1::axpby( 0.5*nu_parallel[i], lambda, 1., lambda,lambda); 
-            }           
-        }
-        dg::blas1::axpby(1.,one,1., logn[i] ,chi); //chi = (1+lnN_e)
-        dg::blas1::axpby(1.,phi[i],p.tau[i], chi); //chi = (tau_e(1+lnN_e)+phi)
-        dg::blas1::pointwiseDot(y[i+2],y[i+2], omega);  
-        dg::blas1::axpby(0.5*p.mu[i], omega,1., chi); //chi = (tau_e(1+lnN_e)+phi + 0.5 mu U^2)
-        if( p.nu_parallel != 0)
-        {
-            Dpar[i] = z[i]*dg::blas2::dot(chi, w3d, lambda); //Z*(tau (1+lnN )+psi) nu_para *(ds^2 N -ds lnB ds N)
-            if( i==0) //only electrons
-            {
-                dg::blas1::axpby(1.,one,1., logn[i] ,chi); //chi = (1+lnN)
-                aligned_ = dg::blas2::dot( chi, w3d, lambda); //(1+lnN)*Delta_s N
-            }
-        }
-        else 
-        {
-            if(i==0) aligned_ = 0;
-            Dpar[i] = 0;
-        }
-
-        //Compute perp dissipation 
-        dg::blas2::gemv( lapperpN, y[i], lambda);
-        dg::blas2::gemv( lapperpN, lambda, omega);//nabla_RZ^4 N_e
-        Dperp[i] = -z[i]* p.nu_perp*dg::blas2::dot(chi, w3d, omega);  
-
-        //Compute parallel dissipative energy for U/////////////////////////////
-        if( p.nu_parallel !=0)
-        {
-            if (p.pollim==1) dsDIR_.set_boundaries( dg::DIR, 0, 0);    
-            if (p.pardiss==0)
-            {
-                dsDIR_.forward(y[i+2], omega); 
-                dsDIR_.forwardT(omega,lambda);
-                dg::blas1::axpby( 0.5*nu_parallel[i+2], lambda, 0.,lambda,lambda); 
-                dsDIR_.backward( y[i+2], omega); 
-                dsDIR_.backwardT(omega,lambda);
-                dg::blas1::axpby( 0.5*nu_parallel[i+2], lambda,  1., lambda,lambda); 
-            }
-            if (p.pardiss==1)
-            {
-/*                    dsDIR_.dsz(y[i+2],omega);                                          //ds^2 U 
-                dg::blas1::axpby( nu_parallel[i+2], omega, 0., lambda,lambda);     //lambda = nu_para*ds^2 U 
-                //gradlnBcorrection
-                dsDIR_(y[i+2], dsy[i+2]);                                               //ds U
-                dg::blas1::pointwiseDot(gradlnB,dsy[i+2], omega);               // ds lnB ds U
-                dg::blas1::axpby(-nu_parallel[i+2], omega, 1., lambda,lambda);     // lambda += nu_para*ds lnB ds N */  
-                dsDIR_.forward( y[i+2], omega); 
-                dsN_.forwardTD(omega,lambda);
-                dg::blas1::axpby( 0.5*nu_parallel[i+2], lambda, 1., lambda,lambda); 
-                dsDIR_.backward( y[i+2], omega); 
-                dsN_.backwardTD(omega,lambda);
-                dg::blas1::axpby( 0.5*nu_parallel[i+2], lambda, 1., lambda,lambda); 
-            }   
-        }
-        dg::blas1::pointwiseDot( npe[i], y[i+2], omega); //N U   
-        if( p.nu_parallel !=0)
-            Dpar[i+2] = z[i]*p.mu[i]*dg::blas2::dot(omega, w3d, lambda);      //Z*N*U nu_para *(ds^2 U -ds lnB ds U)  
-        else
-            Dpar[i+2] = 0;
-
-        //Compute perp dissipation 
-        dg::blas2::gemv( lapperpDIR, y[i+2], lambda);
-        dg::blas2::gemv( lapperpDIR, lambda,chi);//nabla_RZ^4 U
-        Dperp[i+2] = -z[i]*p.mu[i]*p.nu_perp* dg::blas2::dot(omega, w3d, chi);
-
-
-    }
-    //Compute rhs of energy theorem
-    ediff_= Dpar[0]+Dperp[0]+Dpar[1]+Dperp[1]+Dpar[2]+Dperp[2]+Dpar[3]+Dperp[3] + Dres;
-    //end of energy computation
     
     for( unsigned i=0; i<2; i++)
     {

            dg::blas2::gemv( LaplacianM_perpDIR, x[i+2], temp);
            dg::blas2::gemv( LaplacianM_perpDIR, temp, y[i+2]);
            dg::blas1::scal( y[i+2], -p.nu_perp);  //  nu_perp lapl_RZ (lapl_RZ N) 
//             if (p.pardiss==2)
//             {
//             dg::blas2::gemv( dsTdsN, x[i],temp); //lapd is negative
//             dg::blas1::axpby(  -p.nu_parallel ,temp, 1., y[i]);
//             dg::blas2::gemv( dsTdsDIR, x[i+2],temp); //lapd is negative
//             dg::blas1::axpby(  -p.nu_parallel ,temp, 1., y[i+2]);
//             }
        }
        //Resistivity
        dg::blas1::axpby( 1., x[3], -1, x[2], omega); //U_i - U_e

    double z[2]    = {-1.0,1.0};
    double S[2]    = {0.0, 0.0};
    double Tpar[2] = {0.0, 0.0};
    //transform compute n and logn and energies
    for(unsigned i=0; i<2; i++)
    {
