             transferH = dg::pullback( std::get<2>(tp), gX2d);
             dg::blas1::pointwiseDot( volX2d, transferH, transferH);
             avg_eta( transferH, transferH1d, false);
-            dg::blas1::scal( transferH1d, 4*M_PI*M_PI); //
-            dg::blas1::pointwiseDivide( transferH1d, dvdzeta, transferH1d );
+            dg::blas1::scal( transferH1d, 4*M_PI*M_PI*f0); //
+            dg::blas1::pointwiseDivide( transferH1d, dvdpsip, transferH1d );
             map1d.emplace_back( std::get<0>(tp)+"_fsa", transferH1d,
                 std::get<1>(tp)+" (Flux surface average)");
-            dg::blas1::pointwiseDot( transferH1d, dvdzeta, transferH1d );
+            dg::blas1::pointwiseDot( transferH1d, dvdpsip, transferH1d );
             transferH1d = dg::integrate( transferH1d, gX1d);
             map1d.emplace_back( std::get<0>(tp)+"_ifs", transferH1d,
                 std::get<1>(tp)+" (Flux surface integral)");

    };

    /// -------  Elements for fsa on X-point grid ----------------
    if( gp.hasXpoint())
    {
        std::cout << "Generate X-point flux-aligned grid ... \n";

    std::vector<std::tuple<std::string, dg::HVec, std::string> > map1d;
    ///////////TEST CURVILINEAR GRID TO COMPUTE FSA QUANTITIES
    unsigned npsi = 3, Npsi = p.Npsi;//set number of psivalues (NPsi % 8 == 0)
    //Generate list of functions to evaluate
    std::vector< std::tuple<std::string, std::string, dg::geo::CylindricalFunctor >> map{
        {"Psip", "Flux function", mag.psip()},
