     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)
-    double psipmax = dg::blas1::reduce( psipog2d, 0. ,thrust::maximum<double>()); //DEPENDS ON GRID RESOLUTION!!
     //Generate list of functions to evaluate
     std::vector< std::tuple<std::string, std::string, dg::geo::CylindricalFunctor >> map{
         {"Psip", "Flux function", mag.psip()},

    }
    /// --------- More flux labels --------------------------------
    dg::Grid1d grid1d(psipmin<psipmax ? psipmin : psipmax,
            psipmin<psipmax ? psipmax : psipmin, npsi ,Npsi,dg::NEU);
    if( !gp.isToroidal())
    {
        dg::HVec rho = dg::evaluate( dg::cooX1d, grid1d);

        dg::SparseTensor<dg::HVec> metricX = gX2d.metric();
        dg::HVec volX2d = dg::tensor::volume2d( metricX);
        dg::blas1::pointwiseDot( coordsX[0], volX2d, volX2d); //R\sqrt{g}
        const double f0 = (gX2d.x1()-gX2d.x0())/ ( psipmax - psipO);
        dg::HVec dvdpsip;
        avg_eta( volX2d, dvdpsip, false);
        dg::blas1::scal( dvdpsip, 4.*M_PI*M_PI*f0);
        dg::Grid1d gX1d(psipmin<psipmax ? psipmin : psipmax,
            psipmin<psipmax ? psipmax : psipmin, npsi ,Npsi,dg::DIR_NEU); //inner value is always zero
        dg::HVec X_psi_vol = dg::integrate( dvdpsip, gX1d);
        map1d.emplace_back( "dvdpsip", dvdpsip,
            "Derivative of flux volume with respect to flux label psi");
        map1d.emplace_back( "psi_vol", X_psi_vol,
            "Flux volume on X-point grid");

