     dg::DVec temp0( g2d.size()), temp1(temp0);
     dg::DVec w3d = dg::create::weights( g2d);
 
-    //err = nc_put_var_double( ncid, coordsID[0], periodify(X, g2d_periodic).data());
-    //err = nc_put_var_double( ncid, coordsID[1], periodify(Y, g2d_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, 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(), g2d_periodic).data());
-    err = nc_put_var_double( ncid, coord1D[4], g3d.f_x().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, g2d_periodic).data());
-    err = nc_put_var_double( ncid, defID, X.data());
+    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, g2d_periodic).data());
-    err = nc_put_var_double( ncid, volID, X.data());
+    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";
 

double cosine( double x) {return cos(x);}
typedef dg::FieldAligned< orthogonal::GridX3d<dg::DVec> , dg::IDMatrix, dg::DVec> DFA;

thrust::host_vector<double> periodify( const thrust::host_vector<double>& in, const dg::GridX2d& g)
{
    thrust::host_vector<double> out(g.size());
    for( unsigned i=0; i<g.Ny()-1; i++)
    for( unsigned k=0; k<g.n(); k++)
    for( unsigned j=0; j<g.Nx(); j++)
    for( unsigned l=0; l<g.n(); l++)
        out[((i*g.n() + k)*g.Nx() + j)*g.n()+l] = 
            in[((i*g.n() + k)*g.Nx() + j)*g.n()+l];
    for( unsigned i=g.Ny()-1; i<g.Ny(); i++)
    for( unsigned k=0; k<g.n(); k++)
    for( unsigned j=0; j<g.inner_Nx(); j++)
    for( unsigned l=0; l<g.n(); l++)
        out[((i*g.n() + k)*g.Nx() + j)*g.n()+l] = 
            in[((0*g.n() + k)*g.Nx() + j)*g.n()+l];
    //////////////////////////////////
    for( unsigned i=g.Ny()-1; i<g.Ny(); i++)
    for( unsigned k=0; k<g.n(); k++)
    for( unsigned j=g.inner_Nx(); j<g.Nx(); j++)
    for( unsigned l=0; l<g.n(); l++)
        out[((i*g.n() + k)*g.Nx() + j)*g.n()+l] = 
            in[(((i-1)*g.n() + k)*g.Nx() + j)*g.n()+l];
    return out;
}


    int ncid;
    file::NC_Error_Handle err;
    err = nc_create( "testX.nc", NC_NETCDF4|NC_CLOBBER, &ncid);
    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, 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, 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();
