         dg::blas1::pointwiseDot( damping, yp[i+2], yp[i+2]); 
     }
     //add particle source to dtN
-    //dtN_e
-    dg::blas1::axpby(+1.0,profne,-1.0,npe[0],lambda);//lambda = ne0 - ne    
-    dg::blas1::pointwiseDot(source,lambda,omega);//tanhSource on profNe
-    dg::blas1::transform(omega, omega, dg::POSVALUE<double>()); 
-    dg::blas1::axpby(p.omega_source,omega,1.0,yp[0]);
-    //dtN_i
-    dg::blas1::axpby(p.omega_source,omega,1.0,yp[1]);
-    //add FLR correction
-    dg::blas1::pointwiseDot(source,lambda,lambda);//tanhSource on profNe
-    dg::blas1::transform(lambda, lambda, dg::POSVALUE<double>()); 
-    dg::blas2::gemv( lapperpN, lambda, omega); 
-    dg::blas1::axpby(-p.omega_source*0.5*p.tau[1]*p.mu[1],omega,1.0,yp[1]);   
-
+    if (p.omega_source!=0.0)
+    {
+      //dtN_e
+      dg::blas1::axpby(+1.0,profne,-1.0,npe[0],lambda);//lambda = ne0 - ne    
+      dg::blas1::pointwiseDot(source,lambda,omega);//tanhSource on profNe
+      dg::blas1::transform(omega, omega, dg::POSVALUE<double>()); 
+      dg::blas1::axpby(p.omega_source,omega,1.0,yp[0]);
+      //dtN_i
+      dg::blas1::axpby(p.omega_source,omega,1.0,yp[1]);
+      //add FLR correction
+      dg::blas1::pointwiseDot(source,lambda,lambda);//tanhSource on profNe
+      dg::blas1::transform(lambda, lambda, dg::POSVALUE<double>()); 
+      dg::blas2::gemv( lapperpN, lambda, omega); 
+      dg::blas1::axpby(-p.omega_source*0.5*p.tau[1]*p.mu[1],omega,1.0,yp[1]);   
+    }
     t.toc();
     #ifdef MPI_VERSION
         int rank;

    dg::blas1::transfer( dg::evaluate( dg::zero, g), lambda ); 
    dg::blas1::transfer( dg::evaluate( dg::one,  g), one); 
    phi.resize(2);apar.resize(2); phi[0] = phi[1] = apar[0]=apar[1] =  chi;
    curvapar.resize(2); curvapar[0] = curvapar[1] = chi;
    curvphi = curvkappaphi = npe = logn = u = u2 = un =  phi;
    dsphi = dslogn = dsun = dsu2 = poissonn = poissonu =  phi;
    poissonu2 =poissonw = poissonun = poissonlogn =poissonphi  = phi;


    for( unsigned i=0; i<2; i++)
    {
        //Here we use K(Apar) instead of K_nablaB (Apar) (?)
        vecdotnablaDIR(curvX, curvY,  apar[i], dsblnB[i]);     //K(apar) = K(apar)
        dg::blas1::axpby(  1.,  dslnB,p.beta,  dsblnB[i],dsblnB[i]);  // ds ln B + beta K(apar) 
        
        //ExB dynamics
        poissonN( y[i],   phi[i], yp[i]);                           //[N-1,phi]_RZ
        poissonDIR( y[i+2], phi[i], yp[i+2]);                       //[w,phi]_RZ 
        dg::blas1::pointwiseDot( yp[i], binv, yp[i]);               // dtN =1/B [N,phi]_RZ
        dg::blas1::pointwiseDot( yp[2+i], binv, yp[2+i]);           // dtw =1/B [w,phi]_RZ  
        
