 #include "dg/functors.h"
 
 #include "dg/backend/timer.cuh"
+#include "dg/geometry/curvilinear.h"
 //#include "guenther.h"
 #include "solovev.h"
 #include "ribeiro.h"
-#include "curvilinear.h"
-#include "refined_curvilinear.h"
 //#include "ds.h"
 #include "init.h"
 

    }
    //write parameters from file into variables
    dg::geo::solovev::GeomParameters gp(js);
    dg::geo::solovev::Psip psip( gp); 
    std::cout << "Psi min "<<psip(gp.R_0, 0)<<"\n";
    std::cout << "Type psi_0 and psi_1\n";
    double psi_0, psi_1;
    std::cin >> psi_0>> psi_1;
    std::cout << "Type new_n, multiple_x and multiple_y \n";
    double n_ref, multiple_x, multiple_y;
    std::cin >> n_ref>>multiple_x >> multiple_y;
    gp.display( std::cout);
    dg::Timer t;
    //solovev::detail::Fpsi fpsi( gp, -10);
    std::cout << "Constructing ribeiro grid ... \n";
    t.tic();
    dg::geo::solovev::MagneticField c( gp);
    dg::geo::Ribeiro<Psip, PsipR, PsipZ, PsipRR, PsipRZ, PsipZZ>
        ribeiro( c.psip, c.psipR, c.psipZ, c.psipRR, c.psipRZ, c.psipZZ, psi_0, psi_1, gp.R_0, 0., 1);
    //dg::CurvilinearGrid3d<dg::HVec> g3d(ribeiro, n, Nx, Ny,Nz, dg::DIR);
    //dg::CurvilinearGrid2d<dg::HVec> g2d = g3d.perp_grid();
    dg::CurvilinearRefinedGrid3d<dg::HVec> g3d(multiple_x, multiple_y, &ribeiro, n_ref, n, Nx, Ny,Nz, dg::DIR);
    dg::CurvilinearRefinedGrid2d<dg::HVec> g2d = g3d.perp_grid();
    dg::Grid2d g2d_periodic(g2d.x0(), g2d.x1(), g2d.y0(), g2d.y1(), g2d.n(), g2d.Nx(), g2d.Ny()+1); 
    t.toc();
    std::cout << "Construction took "<<t.diff()<<"s"<<std::endl;

    err = nc_def_var( ncid, "volume", NC_DOUBLE, 2, dim3d, &volID);
    err = nc_def_var( ncid, "divB", NC_DOUBLE, 2, dim3d, &divBID);

    thrust::host_vector<double> psi_p = dg::pullback( psip.f(), g2d);
    //g.display();
    err = nc_put_var_double( ncid, onesID, periodify(psi_p, g2d_periodic).data());
    dg::HVec X( g2d.size()), Y(X); //P = dg::pullback( dg::coo3, g);
    for( unsigned i=0; i<g2d.size(); i++)
    {
        X[i] = g2d.map()[0][i];
        Y[i] = g2d.map()[1][i];
    }

    dg::HVec temp0( g2d.size()), temp1(temp0);
