 
     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];


    std::vector<double> c(13);
    //read and store geom data
    std::vector<double> v;
    try{ v = file::read_input( "geometry_params.txt"); }
    catch (toefl::Message& m) {  
        m.display(); 
        for( unsigned i = 0; i<v.size(); i++)
            std::cout << v[i] << " ";
            std::cout << std::endl;
        return -1;}
    for( unsigned i = 1; i<v.size(); i++)
    std::stringstream title;
    //write parameters from file into variables
    const solovev::GeomParameters gp(v);
    gp.display( std::cout);

    Rmin=gp.R_0-gp.a;
    Zmin=-gp.a*gp.elongation;
    Rmax=gp.R_0+gp.a; 
    Zmax=gp.a*gp.elongation;
    std::cout << "The grid parameters" <<"\n";
    std::cout  << Rmin<<"rho_s " << Rmax <<"rho_s " << Zmin <<"rho_s " <<Zmax <<"rho_s " <<"\n";
    std::cout << "Type n, Nx, Ny, Nz\n";
    std::cin >> n>> Nx>>Ny>>Nz;
    
    
    solovev::Field field(gp);
    solovev::Psip psip(gp.R_0,gp.A,gp.c);
    solovev::PsipR psipR(gp.R_0,gp.A,gp.c);

    err = nc_put_var_double( ncid, vecID[1], vecZ.data());
    err = nc_put_var_double( ncid, vecID[2], vecP.data());
    nc_close(ncid);
    std::cout << "Check single field by integrating from 0 to 2pi" << "\n";
    dg::HVec v5(1, 0);
    std::vector<thrust::host_vector<double> > in(3, v5);
    std::vector<thrust::host_vector<double> > out(3, v5);
    in[0][0]=gp.R_0+0.9*gp.a; 
//     in[1][0]=0.9*gp.a*gp.elongation;
    in[1][0]=0.0;

    in[2][0]=0.;
    

    dg::integrateRK4( field, in, out,  2*M_PI, gp.rk4eps);
    
    std::cout <<"Rin =  "<< in[0][0] <<" Zin =  "<<in[1][0] <<" sin  = "<<in[2][0]<<"\n";
    std::cout <<"Rout = "<< out[0][0]<<" Zout = "<<out[1][0]<<" sout = "<<out[2][0]<<"\n";

