     const double Zmin=-p.boxscaleZm*mag.params().a();
     const double Rmax=mag.R0()+p.boxscaleRp*mag.params().a();
     const double Zmax=p.boxscaleZp*mag.params().a();
-    const unsigned FACTOR = 6;
 
     unsigned cx = js["output"]["compression"].get(0u,1).asUInt();
     unsigned cy = js["output"]["compression"].get(1u,1).asUInt();

            dg::assign( transferH, transferH_out_float);

            err = nc_enddef( ncid_out);
            err = nc_put_var_float( ncid_out, vID, append(transferH_out_float, g3d_out_equidistant).data());
            err = nc_redef(ncid_out);
        }
    }
    for( auto record : feltor::generate_cyl2cart( g3d_out_equidistant) )
    {
        int vID;
        err = nc_def_var( ncid_out, std::get<0>(record).data(), NC_FLOAT, 3, &dim_ids[1],
            &vID);
        err = nc_put_att_text( ncid_out, vID, "long_name", std::get<1>(record).size(),
            std::get<1>(record).data());
        dg::assign( std::get<2>(record), transferH_out_float);
        err = nc_enddef( ncid_out);
        err = nc_put_var_float( ncid_out, vID, append(transferH_out_float, g3d_out_equidistant).data());
        err = nc_redef(ncid_out);

    }
    for( auto& record : feltor::diagnostics3d_list)
    {
        std::string name = record.name;

    dg::RealCylindricalGrid3d<double> g3d_in( Rmin,Rmax, Zmin,Zmax, 0, 2*M_PI,
        n_out, Nx_out, Ny_out, Nz_out, p.bcxN, p.bcyN, dg::PER);
    dg::RealCylindricalGrid3d<double> g3d_out( Rmin,Rmax, Zmin,Zmax, 0, 2*M_PI,
        n_out, Nx_out, Ny_out, INTERPOLATE*Nz_out, p.bcxN, p.bcyN, dg::PER);
    dg::RealCylindricalGrid3d<double> g3d_out_equidistant( Rmin,Rmax,
            Zmin,Zmax, 0, 2*M_PI, 1, n_out*Nx_out, n_out*Ny_out,
            INTERPOLATE*Nz_out, p.bcxN, p.bcyN, dg::PER);
    dg::RealCylindricalGrid3d<float> g3d_out_periodic( Rmin,Rmax, Zmin,Zmax, 0,
            2*M_PI+g3d_out.hz(), n_out, Nx_out, Ny_out, INTERPOLATE*Nz_out+1,
            p.bcxN, p.bcyN, dg::PER);
    dg::RealCylindricalGrid3d<float> g3d_out_periodic_equidistant( Rmin,Rmax,
            Zmin,Zmax, 0, 2*M_PI+g3d_out.hz(), 1, n_out*Nx_out, n_out*Ny_out,
            INTERPOLATE*Nz_out+1, p.bcxN, p.bcyN, dg::PER);
    // For field-aligned output
    dg::RealCylindricalGrid3d<double> g3d_out_fieldaligned( Rmin,Rmax,
            Zmin,Zmax, -2*M_PI, 2*M_PI, 1, n_out*Nx_out, n_out*Ny_out, 2*Nz_out,
            p.bcxN, p.bcyN, dg::PER);
    dg::RealGrid2d<double> g2d( Rmin,Rmax,Zmin,Zmax, 1, n_out*Nx_out,
            n_out*Ny_out, p.bcxN, p.bcyN);


    // Construct weights and temporaries
    dg::HVec transferH_in = dg::evaluate(dg::zero,g3d_in);
    dg::HVec transferH_out = dg::evaluate(dg::zero,g3d_out);
    dg::HVec transferH = dg::evaluate(dg::zero,g3d_out_equidistant);
    dg::fHVec transferH_out_float = dg::construct<dg::fHVec>( transferH);
    // Construct for field-aligned output
    dg::HVec transferH_aligned_out = dg::evaluate( dg::zero, g3d_out_fieldaligned);
    dg::HVec RR = dg::evaluate( dg::zero, g3d_out_fieldaligned), ZZ(RR), PP(RR);
    std::array<thrust::host_vector<double>,3> yy0{
        dg::evaluate( dg::cooX2d, g2d),
        dg::evaluate( dg::cooY2d, g2d),
        dg::evaluate( dg::zero, g2d)}, yy1(yy0); //s
    auto bhat = dg::geo::createBHat( mag);
    dg::geo::detail::DSFieldCylindrical cyl_field(bhat);
    double eps = 1e-5;
    double deltaPhi = g3d_out_fieldaligned.hz();
    double phi0 = 0.;
    for( unsigned i=0; i<g2d.size(); i++)
    {
        RR[ Nz_out*g2d.size() + i] = yy0[0][i];
        ZZ[ Nz_out*g2d.size() + i] = yy0[1][i];
        PP[ Nz_out*g2d.size() + i] = phi0;
    }
    for( unsigned k=1; k<Nz_out; k++)
    {
        for( unsigned i=0; i<g2d.size(); i++)
        {
            double phi1 = phi0 + deltaPhi;
            std::array<double,3> coords0{yy0[0][i],yy0[1][i],yy0[2][i]}, coords1;
            dg::integrateERK( "Dormand-Prince-7-4-5", cyl_field, phi0, coords0,
                    phi1, coords1, 0., dg::pid_control,
                    dg::geo::detail::ds_norm, eps, 1e-10, g2d );
            yy1[0][i] = coords1[0], yy1[1][i] = coords1[1], yy1[2][i] = coords1[2];
            //now write into right place in RR ...
            RR[ (Nz_out + k)*g2d.size() + i] = yy1[0][i];
            ZZ[ (Nz_out + k)*g2d.size() + i] = yy1[1][i];
            // Note that phi1 can be different from phi0 + deltaPhi
            PP[ (Nz_out + k)*g2d.size() + i] = phi1;
        }
        std::swap( yy0, yy1);
        phi0 += deltaPhi;
    }
    phi0 = 0;
    yy0 = std::array<thrust::host_vector<double>,3> {
        dg::evaluate( dg::cooX2d, g2d),
        dg::evaluate( dg::cooY2d, g2d),
        dg::evaluate( dg::zero, g2d)};
    for( unsigned k=0; k<Nz_out; k++)
    {
        for( unsigned i=0; i<g2d.size(); i++)
        {
            double phi1 = phi0 - deltaPhi;
            std::array<double,3> coords0{yy0[0][i],yy0[1][i],yy0[2][i]}, coords1;
            dg::integrateERK( "Dormand-Prince-7-4-5", cyl_field, phi0, coords0,
                    phi1, coords1, 0., dg::pid_control,
                    dg::geo::detail::ds_norm, eps, 1e-10, g2d);
            yy1[0][i] = coords1[0], yy1[1][i] = coords1[1], yy1[2][i] = coords1[2];
            //now write into right place in RR ...
            RR[ (Nz_out - k - 1)*g2d.size() + i] = yy1[0][i];
            ZZ[ (Nz_out - k - 1)*g2d.size() + i] = yy1[1][i];
            // Note that phi1 can be different from phi0 - deltaPhi
            PP[ (Nz_out - k - 1)*g2d.size() + i] = phi1 + 2*M_PI;
        }
        std::swap( yy0, yy1);
        phi0 -= deltaPhi;
    }
    dg::IHMatrix big_matrix
        = dg::create::interpolation( RR, ZZ, PP, g3d_in);


    // define 4d dimension
    int dim_ids[4], tvarID;
    err = dg::file::define_dimensions( ncid_out, dim_ids, &tvarID,
            g3d_out_periodic_equidistant, {"time", "z", "y", "x"});
    int dim_idsF[4], tvarIDF;
    err = dg::file::define_dimensions( ncid_out, dim_idsF, &tvarIDF,
            g3d_out_fieldaligned, {"timef", "zf", "yf", "xf"});
    std::map<std::string, int> id4d;

    /////////////////////////////////////////////////////////////////////////
    dg::geo::Fieldaligned<dg::CylindricalGrid3d, dg::IHMatrix, dg::HVec> fieldaligned(
        bhat, g3d_out, dg::NEU, dg::NEU, dg::geo::NoLimiter(),
        //let's take NEU bc because N is not homogeneous
        p.rk4eps, 5, 5);
    dg::IHMatrix interpolate_in_2d = dg::create::interpolation(
            g3d_out_equidistant, g3d_out);


    for( auto& record : feltor::diagnostics3d_static_list)
