anchor.diff:          do q=LL,UU
             ! no vect
             CC(1)=CC(1)+(EE(c1,q+N/2,1)+k1(q+N/2,1)/2)*(EE(c1,q-p(pp)+N/2,1)+k1(q-p(pp)+N/2,1)/2)+ &
  -                (EE(c1,q+N/2,2)+k1(q+N/2,2)/2)*(EE(c1,q-p(pp)+N/2,2)+k1(q-p(pp)+N/2,2)/2);
  +                (EE(c1,q+N/2,2)+k1(q+N/2,2)/2)*(EE(c1,q-p(pp)+N/2,2)+k1(q-p(pp)+N/2,2)/2)
             CC(2)=CC(2)+(EE(c1,q+N/2,2)+k1(q+N/2,2)/2)*(EE(c1,q-p(pp)+N/2,1)+k1(q-p(pp)+N/2,1)/2)- &
  -                (EE(c1,q+N/2,1)+k1(q+N/2,1)/2)*(EE(c1,q-p(pp)+N/2,2)+k1(q-p(pp)+N/2,2)/2);
  +                (EE(c1,q+N/2,1)+k1(q+N/2,1)/2)*(EE(c1,q-p(pp)+N/2,2)+k1(q-p(pp)+N/2,2)/2)
           end do
   
           kn2(:)=Tstep*(vv(c1,pp,:)+kv1(:)/2+SSn(:)/2*Tstep)
  
positive.code: 
        do pp=1,N/2
          call random_number(rdist(:2))
          SSn(:)= rdist(:2)*Source_factor_n(pp)/sqrt(Tstep);
  
          LL= max(p(pp)-N/3,-N/3);
          UU= min(N/3,p(pp)+N/3);
          CC(:)=0.0_wp
  
          do q=LL,UU
            ! no vect
            CC(1)=CC(1)+EE(c1,q+N/2,1)*EE(c1,q-p(pp)+N/2,1)+EE(c1,q+N/2,2)*EE(c1,q-p(pp)+N/2,2);
            CC(2)=CC(2)+EE(c1,q+N/2,2)*EE(c1,q-p(pp)+N/2,1)-EE(c1,q+N/2,1)*EE(c1,q-p(pp)+N/2,2);
          end do
  
          kn1(:)=Tstep*(vv(c1,pp,:))
  
negative.code:       gamal2=(-1)*sqrt(pi/8)* omegae / k(ii) / tetabeam(beamj)**2.0_wp * sign(1.0_wp,k(ii)) * nbeam(beami) / &
                n0*omegaL(ii)*(omegaL(ii)-k(ii)*vbeam(beamj)) / (k(ii)*tetabeam(beamj)) * &
                exp((-1)* omegaL(ii)-k(ii)*vbeam(beamj) / k(ii) / tetabeam(beamj)**2.0_wp/2) !Landau damping due to the beam
  !gamal2=(-1)*sqrt(pi/8)*pow(omegae/k(ii)/tetabeam.at(beamj),2)* sign(1,k(ii)) *nbeam.at(beami)/n0*omegaL(ii)*(omegaL(ii)-k(ii)*vbeam.at(beamj))/(k(ii)*tetabeam.at(beamj))*exp((-1)*pow((omegaL(ii)-k(ii)*vbeam.at(beamj))/k(ii)/tetabeam.at(beamj),2)/2)  //Landau damping due to the beam
        gamal3=(-1)*sqrt(pi)*omegae*omegaL(ii)**2.0_wp / k(ii)**3.0_wp * sign(1.0_wp,k(ii)) *se_cte * &
                (1+ omegaL(ii)**2.0_wp/kappa/ k(ii)*theta_se**2.0_wp)**((-1)*(kappa+1))
  
  