             dg::blas2::gemv( LaplacianM_perp, x[i], temp);
             dg::blas2::gemv( LaplacianM_perp, temp, y[i]);
             dg::blas1::axpby( -p.nu_perp, y[i], 0., y[i]); // - nu_perp lapl_RZ (lapl_RZ (lnN,U)) //factor MISSING!?!
+            
+            dg::blas1::pointwiseDot( lapiris_, temp, temp); //N_i U_i
+            dg::blas1::axpby( -0.1, temp, 1., y[i]); // - nu_perp lapl_RZ (lapl_RZ (lnN,U)) //factor MISSING!?!
+
 //             dg::blas2::gemv( LaplacianM_para, x[i], temp);
 //             dg::blas1::axpby(  -p.nu_parallel, temp, 1., y[i]);
         }

    const container& polarisation( const std::vector<container>& y);

    container chi, omega;
    const container binv, curvR, curvZ, gradlnB;
    const container iris, source, damping;

    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 container w3d, v3d, one;

        w3d_( dg::create::w3d(g)), v3d_(dg::create::v3d(g)),
        temp( g.size()),
//         pupil_( dg::evaluate( solovev::Pupil( gp), g))
        lapiris_( dg::evaluate( solovev::TanhDampingInv(gp ), g)),
        pupil_( dg::evaluate( solovev::TanhDamping(gp ), g))
    {
        LaplacianM_perp = dg::create::laplacianM_perp( g, dg::normed, dg::symmetric);
//         LaplacianM_para = dg::create::laplacianM_parallel( g, dg::PER);
