 
     file::NC_Error_Handle err;
     int ncid, dim_ids[3];
-    err = nc_create( "geometry.nc", NC_CLOBBER, &ncid);
-//         err = nc_create( "geometry.nc", NC_NETCDF4|NC_CLOBBER, &ncid);
+//     err = nc_create( "geometry.nc", NC_CLOBBER, &ncid);
+    err = nc_create( "geometry.nc", NC_NETCDF4|NC_CLOBBER, &ncid);
 
     err = file::define_dimensions( ncid, dim_ids, grid);
     int vecID[3];

    double R_0;
}; 

int main()
{
    double Rmin,Zmin,Rmax,Zmax;
//     double A,R_0,a, elongation;

                std::cout << "Rel Diff = "<<reldiff3 <<"\n";
                dzerrfile1 << pow(2,zz)*Nz <<" " << reldiff << std::endl;
                dzerrfile2 << pow(2,zz)*Nz <<" " << reldiff2 << std::endl;
                
                std::cout <<"----(4) test div(B) != 0 "<<"\n";
                dg::DVec bRongrid = dg::evaluate( fieldR, grid);
                dg::DVec bZongrid = dg::evaluate( fieldZ, grid);
                dg::DVec dRbR(g3d.size());
                dg::DVec dZbZ(g3d.size());                
                dg::DVec invRbR(g3d.size());                
//                    //cut boundaries
//                 dg::blas1::pointwiseDot( pupilongrid, bRongrid,bRongrid); 
//                 dg::blas1::pointwiseDot( pupilongrid, bZongrid, bZongrid); 
                
                dg::DVec divB(g3d.size());                
                dg::blas2::gemv( arakawa.dx(), bRongrid, dRbR);
                dg::blas2::gemv( arakawa.dy(), bZongrid, dZbZ);
                dg::blas1::pointwiseDot( invnormrongrid , bRongrid, invRbR);
                dg::blas1::axpby( 1., dRbR   , 1., dZbZ, divB);
                dg::blas1::axpby( 1./gp.R_0, invRbR , 1., divB);
                double normdivB2= dg::blas2::dot( w3d, divB); 
                std::cout << "divB = "<<sqrt( normdivB2)<<"\n";
                
             }
            dzerrfile1.close();
            dzerrfile2.close();
