     err = nc_create( "test.nc", NC_NETCDF4|NC_CLOBBER, &ncid);
     int dim3d[3];
     err = file::define_dimensions(  ncid, dim3d, g.grid());
-    int coordsID[2], onesID;
-    err = nc_def_var( ncid, "r_xy", NC_DOUBLE, 3, dim3d, &coordsID[0]);
-    err = nc_def_var( ncid, "z_xy", NC_DOUBLE, 3, dim3d, &coordsID[1]);
+    int coordsID[3], onesID;
+    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, "z_XYP", NC_DOUBLE, 3, dim3d, &coordsID[2]);
     err = nc_def_var( ncid, "psi", NC_DOUBLE, 3, dim3d, &onesID);
 
-    thrust::host_vector<double> ones = dg::pullback( psip, g);
+    thrust::host_vector<double> psi_p = dg::pullback( psip, g);
     g.grid().display();
-    err = nc_put_var_double( ncid, onesID, ones.data());
-    err = nc_put_var_double( ncid, coordsID[0], g.r().data());
-    err = nc_put_var_double( ncid, coordsID[1], g.z().data());
+    err = nc_put_var_double( ncid, onesID, psi_p.data());
+    dg::HVec X( g.grid().size()), Y(X), P = dg::pullback( dg::coo3, g);
+    for( unsigned i=0; i<g.grid().size(); i++)
+    {
+        X[i] = g.r()[i]*cos(P[i]);
+        Y[i] = g.r()[i]*sin(P[i]);
+    }
+
+
+    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[2], g.z().data());
     err = nc_close( ncid);
 
     std::cout << "Construction successful!\n";
 
     //compute error in volume element
+    dg::HVec ones = dg::evaluate( dg::one, g.grid());
     dg::HVec temp0( g.grid().size()), temp1(temp0);
     dg::HVec w3d = dg::create::weights( g.grid());
     dg::blas1::pointwiseDot( g.g_xx(), g.g_yy(), temp0);

    std::cout << "volumeXYP is "<< volume<<std::endl;
    std::cout << "volumeRZP is "<< volumeRZP<<std::endl;
    std::cout << "relative difference in volume is "<<fabs(volumeRZP - volume)/volume<<std::endl;


    return 0;

    dg::blas1::axpby( 1., temp0, -1., by, temp1);
    double error= dg::blas2::dot( temp1, w3d, temp1);
    std::cout << "Rel Error of g.g_xx() is "<<sqrt(error/dg::blas2::dot( by, w3d, by))<<"\n";
    double volume = dg::blas1::dot( g.vol(), w3d);

    std::cout << "TEST VOLUME IS:\n";
