   err_out = nc_inq_varid(ncid_out, "y", &var_id);
   err_out = nc_put_var_double(ncid_out, var_id, transfer_y);
   //timestepping
+  bool reset_flag[3] = {true, true, true};
+  int m[3] = {};
   int add_to_spatial_idx = 0;
   for (unsigned i = 0; i <= p.maxout; i++)
   { std::cout << i << "\n";   // tbr!

//dg::Timer t;
//t.tic();
//t.toc();
// boost::timer t;
// double duration;
//
// t.restart();
// duration = t.elapsed();

struct Heaviside2d
{ Heaviside2d( double sigma):sigma2_(sigma*sigma), x_(0), y_(0) {}

  //calculation variables per species
  double mass_[num_species] = {}, cn[num_species] = {};
  double posX = 0, posY =0 ;
  double NposX = 0, NposY =0 ; //
  double posX_init[num_species] = {}, posY_init[num_species] = {};
  double NposX_init[num_species] = {}, NposY_init[num_species] = {}; //
  double posX_old[num_species] ={}, posY_old[num_species] = {};
  double NposX_old[num_species] ={}, NposY_old[num_species] = {};
  double posX_max = 0, posY_max = 0;
  double posX_max_old[num_species] = {}, posY_max_old[num_species] = {};
  double posX_max_hs = 0, posY_max_hs = 0;
  double velX[num_species] = {}, velY[num_species] = {};
  double NvelX[num_species] = {}, NvelY[num_species] = {};
  double velX_old[num_species] = {}, velY_old[num_species] = {};
  double NvelX_old[num_species] = {}, NvelY_old[num_species] = {};
  double velX_max = 0, velY_max = 0;
  double velCOM = 0;
  double accX = 0, accY = 0;
