 
     container chi, omega,lambda;//1d container
 
-    container binv,  curvX, curvY, curvKappaX, curvKappaY,divCurvKappa, gradlnB;
+    container binv,  curvX, curvY, curvKappaX, curvKappaY,divCurvKappa, dslnB;
     container source, damping, one;
     container profne, profNi;
     container w3d, v3d;
     
-    std::vector<container> phi,apar,curvapar, u, u2, npe, logn,un,curvphi,curvkappaphi, dsphi, dslogn,dsun,dsu2; //2d container
+    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
 

        }   

        //compute omega = NU
        dg::blas1::pointwiseDot( npe[i], u[i], omega); //N U   
        Dpar[i+2] = z[i]*p.mu[i]*dg::blas2::dot(omega, w3d, lambda);      //Z*N*U nu_para *(ds^2 U ) or Z*N*U nu_para *(ds^2 w )

        //Compute perp dissipation  for U
        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);

    }

{   
    /*  y[0] := N_e - 1
        y[1] := N_i - 1
        y[2] := w_e =  U_e + beta/mu_e Apar_e
        y[3] := w_i =  U_i + beta/mu_i Apar_i
    */
    
    dg::Timer t;
