             std::cerr << "Continue with next file\n";
             continue;
         }
-        err = nc_inq_unlimdim( ncid, &timeID); //Attention: Finds first unlimited dim, which hopefully is time and not energy_time
+        err = nc_inq_unlimdim( ncid, &timeID);
         err = nc_inq_dimlen( ncid, timeID, &steps);
         //steps = 3;
         for( unsigned i=0; i<steps; i++)//timestepping

    err = nc_inq_attlen( ncid_in, NC_GLOBAL, "inputfile", &length);
    std::string inputfile(length, 'x');
    err = nc_get_att_text( ncid_in, NC_GLOBAL, "inputfile", &inputfile[0]);
    err = nc_inq_attlen( ncid_in, NC_GLOBAL, "geomfile", &length);
    std::string geomfile(length, 'x');
    err = nc_get_att_text( ncid_in, NC_GLOBAL, "geomfile", &geomfile[0]);
    err = nc_close( ncid_in);
    Json::Value js,gs;
    dg::file::string2Json(inputfile, js, dg::file::comments::are_forbidden);
    dg::file::string2Json(geomfile, gs, dg::file::comments::are_forbidden);
    //we only need some parameters from p, not all
    const feltor::Parameters p(js, dg::file::error::is_warning);
    const dg::geo::solovev::Parameters gp(gs);
    std::cout << js<<std::endl;
    std::cout << gs<<std::endl;
    std::vector<std::string> names_input{
        "electrons", "ions", "Ue", "Ui", "potential", "aparallel"
    };

    //std::cin >> npsi >> Npsi >> Neta;
    std::cout << "You typed "<<npsi<<" x "<<Npsi<<" x "<<Neta<<"\n";
    std::cout << "Generate X-point flux-aligned grid!\n";
    double R_X = mag.R0()-1.1*mag.params().triangularity()*mag.params().a();
    double Z_X = -1.1*mag.params().elongation()*mag.params().a();
    dg::geo::findXpoint( mag.get_psip(), R_X, Z_X);
    dg::geo::CylindricalSymmTensorLvl1 monitor_chi =
        dg::geo::make_Xconst_monitor( mag.get_psip(), R_X, Z_X) ;
    double R_O = mag.R0(), Z_O = 0;
    dg::geo::findOpoint( mag.get_psip(), R_O, Z_O);
    double psipO = mag.psip()(R_O, Z_O);

    dg::geo::SeparatrixOrthogonal generator(mag.get_psip(), monitor_chi, psipO,
            R_X, Z_X, mag.R0(), 0, 0, false);
    double fx_0 = 1./8.;
    double psipmax = dg::blas1::reduce( psipog2d, 0. ,thrust::maximum<double>()); //DEPENDS ON GRID RESOLUTION!!
    std::cout << "psi max is            "<<psipmax<<"\n";
