     double psipmax = dg::blas1::reduce( psipog2d, 0. ,thrust::maximum<double>()); //DEPENDS ON GRID RESOLUTION!!
     double volumeXGrid;
     /// -------  Elements for fsa of curvature operators ----------------
-    // This is a numerical test because the fsa of the curvature operators exactly vanishes (at least for I = cte)
+    // This is a numerical test because the fsa of the curvature operators exactly vanishes (at least for I = cte and almost for everything else)!!
     dg::geo::CylindricalVectorLvl0 bhat_ = dg::geo::createBHat( mag);
     dg::geo::CylindricalVectorLvl0 curvB_ = dg::geo::createTrueCurvatureNablaB( mag);
     dg::geo::CylindricalVectorLvl0 curvK_ = dg::geo::createTrueCurvatureKappa( mag);

    err = nc_redef(ncid);
    //////////////////////////////Finalize////////////////////////////////////
    err = nc_close(ncid);
    std::cout << "TEST ACCURACY OF CURVATURES (values must be close to 0!)\n";
    std::array<dg::HVec, 3> bhat, curvB, curvK;
    dg::pushForward( bhat_.x(), bhat_.y(), bhat_.z(),
            bhat[0], bhat[1], bhat[2], grid3d);

        double fx_0 = 1./8.;
        psipmax = -fx_0/(1.-fx_0)*psipmin;
        std::cout << "psi 1 is          "<<psipmax<<"\n";
        dg::geo::CurvilinearGridX2d gX2d( generator, fx_0, 0., npsi, Npsi, 640, dg::DIR, dg::NEU);
        std::cout << "DONE! Generate X-point flux-aligned grid!\n";
        dg::Average<dg::HVec > avg_eta( gX2d.grid(), dg::coo2d::y);
        std::vector<dg::HVec> coordsX = gX2d.map();
