         Dperp[i] = -z[i]* p.nu_perp*dg::blas2::dot(chi, w3d, omega);  
 
         dg::blas2::symv(dsDIR_, y[i+2],lambda);
-        dg::blas1::axpby( nu_parallel[i+2], lambda,  0., lambda,lambda); 
+        dg::blas1::axpby( nu_parallel[i+2], lambda,  0., lambda); //nu_para *(ds^2 w )
 
         //compute omega = NU
-        Dpar[i+2] = z[i]*p.mu[i]*dg::blas2::dot(un[i], w3d, lambda);      //Z*N*U nu_para *(ds^2 U ) or Z*N*U nu_para *(ds^2 w )
+        Dpar[i+2] = z[i]*p.mu[i]*dg::blas2::dot(un[i], w3d, lambda); // Z*N*U nu_para *(ds^2 w )
 
         //Compute perp dissipation  for U
         dg::blas2::gemv( lapperpDIR, y[i+2], lambda);

    container w3d, v3d;
    
    std::vector<container> phi,apar,dsblnB, u, u2, npe, logn,un,curvphi,curvkappaphi, dsphi, dslogn,dsun,dsu2; //2d container
    std::vector<container> poissonn,poissonu,poissonw,poissonun,poissonlogn,poissonphi,poissonu2; //2d container
    std::vector<container> dsy, curvy,curvkappay;  //4d container

    //matrices and solvers

    for(unsigned i=0; i<2; i++)
    {
        //compute em Poisson bracket of ds_b
        poissonN( logn[i], apar[i],poissonlogn[i]);                         // -[Apar,logN]_RZ
        poissonDIR( u2[i],   apar[i],poissonu2[i]);                         // -[Apar,U^2]_RZ 
        if (p.bc==dg::NEU)
        {
            poissonN(  y[i],  apar[i], poissonn[i]);                           // -[Apar,N]_RZ 
            poissonDIR( u[i], apar[i], poissonu[i]);                           // -[Apar,U]_RZ
            dg::blas1::pointwiseDot( 1.0, poissonu[i], npe[i], 1.0, poissonn[i],  u[i], 0.0, poissonun[i]); //- N[Apar,U]_RZ  - U[Apar,N]_RZ
        }
        if (p.bc==dg::DIR)
        {
            poissonDIR(un[i], apar[i], poissonun[i]);                         // -[Apar,U N]_RZ 
        }
        //multiply em Poisson bracket by 1/B
        dg::blas1::pointwiseDot(poissonlogn[i], binv, poissonlogn[i]);  //-1/B [Apar,logN]_RZ
        dg::blas1::pointwiseDot(poissonun[i]  , binv, poissonun[i]);    //-1/B [Apar,UN]_RZ
        dg::blas1::pointwiseDot(poissonu2[i]  , binv, poissonu2[i]);    //-1/B [Apar,U^2]_RZ

        //Parallel dynamics (compute ds UN, ds U^2, ds psi, ds ln N) 
        if (p.bc==dg::NEU)
        {
            dsN_.centered(y[i],   chi);                                             // ds N
            dsDIR_.centered(u[i], omega);                                           // ds u
            dg::blas1::pointwiseDot(1.0, chi, u[i], 1.0, omega, npe[i], 0.0, dsun[i]);  // ds UN = N ds U +   U ds N 
            dsN_.centered(logn[i], dslogn[i]);                                            // ds log N
        }
        if (p.bc==dg::DIR)
        {
            dsDIR_.centered(logn[i], dslogn[i]);                                           // ds log N
            dsDIR_.centered(un[i], dsun[i]);                                               // ds U N
        }
        dsDIR_.centered(u2[i],dsu2[i]);                                                  // ds U^2
        dsDIR_.centered(phi[i], dsphi[i]);                                               // ds psi

        //Add terms  ds^b f = ds f-beta/B [Apar,f] terms (Note:  [Apar,phi] term is in ExB contribution  )
        dg::blas1::pointwiseDot(1., un[i], dsblnB[i], 1., yp[i]);                    // dtN += U N ds^b ln B
        dg::blas1::axpbypgz( -1., dsun[i],-p.beta, poissonun[i] ,1., yp[i]);         // dtN += - ds^b U N 
        dg::blas1::axpbypgz( -p.tau[i]/p.mu[i], dslogn[i], -p.tau[i]/p.mu[i]*p.beta, poissonlogn[i], 1., yp[2+i]); // dtw += -tau/hat(mu)*ds^b lnN 
        dg::blas1::axpby( -1./p.mu[i], dsphi[i], 1., yp[2+i]); // dtw += - 1/(hat(mu))  *ds psi   
        dg::blas1::axpbypgz( -0.5, dsu2[i], -0.5*p.beta, poissonu2[i], 1., yp[2+i]); // dtw += - 0.5 ds U^2

        //Parallel dissipation
        double nu_parallel[] = {p.nu_parallel, p.nu_parallel, p.nu_parallel, p.nu_parallel};
