     double Uphi = 0.5*p.mu[1]*dg::blas2::dot( expy[1], w3d, omega); 
     dg::blas1::pointwiseDot( y[2], y[2], omega); //U_e^2
     double Upare = -0.5*p.mu[0]*dg::blas2::dot( expy[0], w3d, omega); 
-    dg::blas1::pointwiseDot( y[3], y[3], omega); //U_i^2
+    dg::blas1::pointwiseDot(y[3], y[3], omega); //U_i^2
     double Upari =  0.5*p.mu[1]*dg::blas2::dot( expy[1], w3d, omega); 
     energy_ = Ue + Ui + Uphi + Upare + Upari;
 
     for( unsigned i=0; i<2; i++)
     {
-
-        arakawa( y[i], phi[i], yp[i]);                      
-        arakawa( y[i+2], phi[i], yp[i+2]);  
-#if def APAR
-        arakawa( apar, y[i], arakAN[i]); // [A_parallel,N]_RZ
-        arakawa( apar, u[i], arakAU[i]); // [A_parallel,U]_RZ
-#endif
-        dg::blas1::pointwiseDot( yp[i], binv, yp[i]);                        // dtlnN =1/B [phi,lnN]_RZ
-        dg::blas1::pointwiseDot( yp[2+i], binv, yp[2+i]);                    // dtU =1/B [phi,U]_RZ
-    
+        //Compute RZ poisson  brackets
+        arakawa( y[i], phi[i], yp[i]);  //[lnN,phi]_RZ
+        arakawa( y[i+2], phi[i], yp[i+2]);//[U,phi]_RZ  
+        dg::blas1::pointwiseDot( yp[i], binv, yp[i]);                        // dtlnN =1/B [lnN,phi]_RZ
+        dg::blas1::pointwiseDot( yp[2+i], binv, yp[2+i]);                    // dtU =1/B [U,phi]_RZ  
+
+
         //compute parallel derivatives
         dz(y[i], dzy[i]);
         dz(phi[i], dzphi[i]);
-        dz(y[2+i], dzy[2+i]);
-    
-        //add A_parallel terms to the parallel derivaitives
-//         dg::axpby(-p.beta,araAparaN ,1.,dzy[i]);
-//         dg::axpby(-p.beta,araAparaphi ,1.,dzphi[i]);
-//         dg::axpby(-p.beta,araAparaU,1.,dzy[2+i]);
+        dz(y[i+2], dzy[2+i]);
+
         //parallel advection terms
-        dg::blas1::pointwiseDot(y[2+i], dzy[i], omega); //Udz lnN 
+        dg::blas1::pointwiseDot(y[i+2], dzy[i], omega); //Udz lnN 
         dg::blas1::axpby( -1., omega, 1., yp[i]);                            //dtlnN = dtlnN - Udz lnN
         dg::blas1::axpby( -1., dzy[2+i], 1., yp[i]);                         //dtlnN = dtlnN - dz U
-        dg::blas1::pointwiseDot(y[2+i], gradlnB, omega);                     //dtlnN = dtlnN + U dz ln B
-        dg::blas1::axpby( 1., omega, 1., yp[i]); 
-        dg::blas1::pointwiseDot(y[2+i], dzy[2+i], omega);                    //dtU = dtU - U dz U
-        dg::blas1::axpby( -1., omega, 1., yp[2+i]);
+        dg::blas1::pointwiseDot(y[i+2], gradlnB, omega);                     
+        dg::blas1::axpby( 1., omega, 1., yp[i]);                            //dtlnN = dtlnN + U dz ln B
+        dg::blas1::pointwiseDot(y[i+2], dzy[2+i], omega);                    
+        dg::blas1::axpby( -1., omega, 1., yp[2+i]);                         //dtw = dtU - U dz U
 
         //parallel force terms
-        dg::blas1::axpby( -p.tau[i]/p.mu[i]/p.eps_hat, dzy[i], 1., yp[2+i]); //dtU = dtU - tau/(hat(mu))*dz lnN
-        dg::blas1::axpby( -1./p.mu[i]/p.eps_hat, dzphi[i], 1., yp[2+i]);     //dtU = dtU - 1/(hat(mu))*dz phi
+        dg::blas1::axpby( -p.tau[i]/p.mu[i]/p.eps_hat, dzy[i], 1., yp[2+i]); //dtw = dtU - tau/(hat(mu))*dz lnN
+        dg::blas1::axpby( -1./p.mu[i]/p.eps_hat, dzphi[i], 1., yp[2+i]);     //dtw = dtU - 1/(hat(mu))*dz phi
 
         //curvature terms
-        curve( y[i], curvy[i]); //K(N)
-        curve( y[2+i], curvy[2+i]);//K(U)
-        curve( phi[i], curvphi[i]);//K(phi)
-
-        dg::blas1::pointwiseDot( y[2+i], curvy[2+i], omega); //U K(U)
-        dg::blas1::pointwiseDot( y[2+i], omega, chi); //U^2 K(U)
+        curve( y[i], curvy[i]);     //K(N)
+        curve( y[i+2], curvy[2+i]); //K(U)
+        curve( phi[i], curvphi[i]); //K(phi)
+
+        dg::blas1::pointwiseDot(y[i+2], curvy[2+i], omega); //U K(U)
+        dg::blas1::pointwiseDot( y[i+2], omega, chi); //U^2 K(U)
         dg::blas1::axpby( -p.mu[i]*p.eps_hat, omega, 1., yp[i]);             //dtlnN = dtlnN - (hat(mu)) U K(U)
-        dg::blas1::axpby( -0.5*p.mu[i]*p.eps_hat, chi, 1., yp[2+i]);         //dtU = dtU - 0.5 (hat(mu)) U^2 K(U)
-
-        dg::blas1::pointwiseDot( y[2+i], curvy[i], omega); //U K(ln N)
-        dg::blas1::pointwiseDot( y[2+i], omega, chi); //U^2K(ln N)
-        dg::blas1::axpby( -p.tau[i], omega, 1., yp[2+i]);                    //dtU = dtU - tau U K(lnN)
+        dg::blas1::axpby( -0.5*p.mu[i]*p.eps_hat, chi, 1., yp[2+i]);         //dtw = dtU - 0.5 (hat(mu)) U^2 K(U)
+
+        dg::blas1::pointwiseDot(y[i+2], curvy[i], omega); //U K(ln N)
+        dg::blas1::pointwiseDot( y[i+2], omega, chi); //U^2K(ln N)
+        dg::blas1::axpby( -p.tau[i], omega, 1., yp[2+i]);                    //dtw = dtU - tau U K(lnN)
         dg::blas1::axpby( -0.5*p.mu[i]*p.eps_hat, chi, 1., yp[i]);           //dtlnN = dtlnN -0.5mu U^2K(lnN)
 
         dg::blas1::axpby( -p.tau[i], curvy[i], 1., yp[i]);                   //dtlnN = dtlnN - tau K(lnN)
-        dg::blas1::axpby( -2.*p.tau[i], curvy[2+i], 1., yp[2+i]);            //dtU = dtU - 2 tau K(U)
+        dg::blas1::axpby( -2.*p.tau[i], curvy[2+i], 1., yp[2+i]);            //dtw = dtU - 2 tau K(U)
         dg::blas1::axpby( -1., curvphi[i], 1., yp[i]);                       //dtlnN= dtlnN - K(psi)
 
-        dg::blas1::pointwiseDot( y[2+i], curvphi[i], omega);  //U K(phi)
-        dg::blas1::axpby( -0.5, omega, 1., yp[2+i]);                         //dtU = dtU -0.5 UK(psi)
-
-
-    }
-    //add parallel resistivity
-//     dg::blas1::pointwiseDot( expy[0], y[2], omega);
-//     dg::blas1::pointwiseDot( expy[1], y[3], chi);
-//     dg::blas1::axpby( -1., omega, 1., chi); //-N_eU_e + N_iU_i
-//     dg::blas1::pointwiseDivide( chi, expy[0], omega);//J_par/N_e
-//     dg::blas1::pointwiseDivide( chi, expy[1], chi); //J_par/N_i
-// 
-//     dg::blas1::axpby( -p.c/p.mu[0]/p.eps_hat, omega, 1., yp[2]);
-//     dg::blas1::axpby( -p.c/p.mu[1]/p.eps_hat, chi, 1., yp[3]);
-
+        dg::blas1::pointwiseDot( y[i+2], curvphi[i], omega);  //U K(phi)
+        dg::blas1::axpby( -0.5, omega, 1., yp[2+i]);                         //dtw = dtU -0.5 U K(psi)
+    }
+        
     //add parallel diffusion with naive implementation
     for( unsigned i=0; i<4; i++)
     {

template< class container>
Feltor< container>::Feltor( const dg::Grid3d<value_type>& g, Parameters p, solovev::GeomParameters gp): 
    chi( g.size(), 0.), omega(chi),
#if def APAR
    rho( grid.size(), 0.),    
#endif
    binv( dg::evaluate(solovev::Field(gp) , g) ),
    curvR( dg::evaluate( solovev::CurvatureR(gp), g)),
    curvZ( dg::evaluate(solovev::CurvatureZ(gp), g)),

    const container& compute_vesqr( container& potential);
    //extrapolates and solves for phi[1], then adds square velocity ( omega)
    const container& compute_psi( container& potential);
    const container& polarisation( const std::vector<container>& y); //solves polarisation equation

    container chi, omega;

    const container binv, curvR, curvZ, gradlnB;
    const container iris, source, damping;
    const container w3d, v3d, one;
    std::vector<container> phi, curvphi, dzphi, expy;
    std::vector<container> dzy, curvy; 

    //matrices and solvers
    Matrix A; 
    dg::DZ<container> dz;
    dg::ArakawaX< dg::DMatrix, container>    arakawa; 
    //dg::Polarisation2dX< thrust::host_vector<value_type> > pol; //note the host vector
    dg::Polarisation< dg::DMatrix, container, container > pol; //note the host vector
    dg::Invert<container> invert_pol;

    const Parameters p;
    const solovev::GeomParameters gp;

