 
 
 template<class DS, class M, class V, class P>
-void Feltor<DS, M, V, P>::add_parallel_dynamics( std::vector<V>& y, std::vector<V>& yp)
-{
+double Feltor<DS, M, V, P>::add_parallel_dynamics( std::vector<V>& y, std::vector<V>& yp)
+{
+    //parallel and perpendicular energy dissipation
+    double Dpar[4] = {0.0, 0.0,0.0,0.0};
+    double Dperp[4] = {0.0, 0.0,0.0,0.0};
+    if (p.pollim==1) dsN_.set_boundaries( p.bc, 0, 0);  //ds N  on limiter
+    if (p.pollim==1) dsDIR_.set_boundaries( dg::DIR, 0, 0); //ds psi on limiter
     for(unsigned i=0; i<2; i++)
     {
         //Parallel dynamics
-        if (p.pollim==1) dsN_.set_boundaries( p.bc, 0, 0);  //ds N  on limiter
         dsN_(y[i], chi);   
         dg::blas1::pointwiseDot(y[i+2], chi, omega);     // U ds N
-        if (p.pollim==1) dsDIR_.set_boundaries( dg::DIR, 0, 0);  //ds U  on limiter
         dsDIR_(y[i+2], chi);  
         dg::blas1::pointwiseDot(npe[i], chi,chi);     // N ds U
         dg::blas1::axpby(1.0,chi,1.0,omega,chi);

            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

    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)
    
    for( unsigned i=0; i<2; i++)
    {
