     int ncid;
     file::NC_Error_Handle err;
     err = nc_create( "testX.nc", NC_NETCDF4|NC_CLOBBER, &ncid);
-    int dim3d[2], dim1d[1];
-    //err = file::define_dimensions(  ncid, dim3d, g2d_periodic.grid());
-    err = file::define_dimensions(  ncid, dim3d, g2d.grid());
+    int dim3d[3], dim1d[1];
+    err = file::define_dimensions(  ncid, dim3d, g3d_periodic.grid());
+    //err = file::define_dimensions(  ncid, dim3d, g2d.grid());
     err = file::define_dimension(  ncid, "i", dim1d, g1d);
     int coordsID[2], onesID, defID, volID, divBID;
     int coord1D[5];
-    err = nc_def_var( ncid, "x_XYP", NC_DOUBLE, 2, dim3d, &coordsID[0]);
-    err = nc_def_var( ncid, "y_XYP", NC_DOUBLE, 2, dim3d, &coordsID[1]);
+    err = nc_def_var( ncid, "x_XYP", NC_DOUBLE, 3, dim3d, &coordsID[0]);
+    err = nc_def_var( ncid, "y_XYP", NC_DOUBLE, 3, dim3d, &coordsID[1]);
     err = nc_def_var( ncid, "x_left", NC_DOUBLE, 1, dim1d, &coord1D[0]);
     err = nc_def_var( ncid, "y_left", NC_DOUBLE, 1, dim1d, &coord1D[1]);
     err = nc_def_var( ncid, "x_right", NC_DOUBLE, 1, dim1d, &coord1D[2]);
     err = nc_def_var( ncid, "y_right", NC_DOUBLE, 1, dim1d, &coord1D[3]);
     err = nc_def_var( ncid, "f_x", NC_DOUBLE, 1, dim1d, &coord1D[4]);
     //err = nc_def_var( ncid, "z_XYP", NC_DOUBLE, 3, dim3d, &coordsID[2]);
-    err = nc_def_var( ncid, "psi", NC_DOUBLE, 2, dim3d, &onesID);
-    err = nc_def_var( ncid, "deformation", NC_DOUBLE, 2, dim3d, &defID);
-    err = nc_def_var( ncid, "volume", NC_DOUBLE, 2, dim3d, &volID);
-    err = nc_def_var( ncid, "divB", NC_DOUBLE, 2, dim3d, &divBID);
+    err = nc_def_var( ncid, "psi", NC_DOUBLE, 3, dim3d, &onesID);
+    err = nc_def_var( ncid, "deformation", NC_DOUBLE, 3, dim3d, &defID);
+    err = nc_def_var( ncid, "volume", NC_DOUBLE, 3, dim3d, &volID);
+    err = nc_def_var( ncid, "divB", NC_DOUBLE, 3, dim3d, &divBID);
 
     thrust::host_vector<double> psi_p = dg::pullback( psip, g2d);
     g2d.display();

    orthogonal::GridX3d<dg::DVec> g3d(gp, psi_0, fx_0, fy_0, n, Nx, Ny,Nz, dg::DIR, dg::NEU);
    orthogonal::GridX2d<dg::DVec> g2d = g3d.perp_grid();
    t.toc();
    //dg::GridX2d g2d_periodic(g2d.x0(), g2d.x1(), g2d.y0(), g2d.y1(), g2d.fx(), g2d.fy(), g2d.n(), g2d.Nx(), g2d.Ny()+1); 
    std::cout << "Construction took "<<t.diff()<<"s"<<std::endl;
    dg::Grid1d<double> g1d( g2d.x0(), g2d.x1(), g2d.n(), g2d.Nx());
    dg::HVec x_left = dg::evaluate( sine, g1d), x_right(x_left);

    dg::DVec temp0( g2d.size()), temp1(temp0);
    dg::DVec w3d = dg::create::weights( g2d);

    err = nc_put_var_double( ncid, coordsID[0], periodify(X, g3d_periodic).data());
    err = nc_put_var_double( ncid, coordsID[1], periodify(Y, g3d_periodic).data());
    //err = nc_put_var_double( ncid, coordsID[0], X.data());
    //err = nc_put_var_double( ncid, coordsID[1], Y.data());
    //err = nc_put_var_double( ncid, coord1D[0], g3d.rx0().data());
    //err = nc_put_var_double( ncid, coord1D[1], g3d.zx0().data());
    //err = nc_put_var_double( ncid, coord1D[2], g3d.rx1().data());
    //err = nc_put_var_double( ncid, coord1D[3], g3d.zx1().data());
    err = nc_put_var_double( ncid, coord1D[4], periodify(g3d.f_x(), g3d_periodic).data());
    //err = nc_put_var_double( ncid, coord1D[4], g3d.f_x().data());
    //err = nc_put_var_double( ncid, coordsID[2], g.z().data());

    dg::blas1::pointwiseDivide( g2d.g_yy(), g2d.g_xx(), temp0);
    dg::blas1::axpby( 1., ones, -1., temp0, temp0);
    X=temp0;
    err = nc_put_var_double( ncid, defID, periodify(X, g3d_periodic).data());
    //err = nc_put_var_double( ncid, defID, X.data());
    X = g2d.vol();
    err = nc_put_var_double( ncid, volID, periodify(X, g3d_periodic).data());
    //err = nc_put_var_double( ncid, volID, X.data());

    std::cout << "Construction successful!\n";

