       err_out = nc_put_vara_double(species_wgrp_id[j], species_wspatial_id[7 + add_to_spatial_idx], start0d, count0d, &posY_max);
       velX_max = (posX_max - posX_max_old[j]) /deltaT;
       velY_max = (posY_max - posY_max_old[j]) /deltaT;
-      if (i>0)
+      if (m[j]>0)
       { posX_max_old[j] = posX_max;
         posY_max_old[j] = posY_max;
       }
       err_out = nc_put_vara_double(species_wgrp_id[j], species_wspatial_id[8 + add_to_spatial_idx], start0d, count0d, &velX_max);
       err_out = nc_put_vara_double(species_wgrp_id[j], species_wspatial_id[9 + add_to_spatial_idx], start0d, count0d, &velY_max);
       //compactness
-      if (i==0)
+      if (m[j]==0)
       { heavi.set_origin( posX_max_hs, posY_max_hs );
         heavy = dg::evaluate( heavi, g2d);
         normalize = dg::blas2::dot( heavy, w2d, ntilde[j]);

  };
  std::string err_time[num_err_time] =
  { "Se", "Si", "Sz",
    "Uphii", "Uphiz", "dcn"
  };
  std::string err_etime[num_err_etime] =
  { "energy", "mass", "dissipation",

  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!
