diff --git a/UWLCM_plotters/include/Plotter3d.hpp b/UWLCM_plotters/include/Plotter3d.hpp
index 9d80fed..4f99ed8 100644
--- a/UWLCM_plotters/include/Plotter3d.hpp
+++ b/UWLCM_plotters/include/Plotter3d.hpp
@@ -18,7 +18,7 @@ class Plotter_t<3> : public PlotterCommon
   using parent_t = PlotterCommon;
   hsize_t n[3];
   enum {x, y, z};
-  arr_t tmp, tmp_srfc;
+  arr_t tmp, tmp_srfc, tmp_ref;
   blitz::Range yrange;
 
   public:
@@ -36,15 +36,17 @@ class Plotter_t<3> : public PlotterCommon
       cnt[3] = { n[x],  n[y],  srfc ? 1 : n[z] }, 
       off[3] = { 0,     0,     0    };
     this->h5s.selectHyperslab( H5S_SELECT_SET, cnt, off);
+
+    arr_t &arr(tmp.extent(0) == n[x] ? tmp : tmp_ref); // crude check if it is refined or normal data; assume surface data is never refined
   
     hsize_t ext[3] = {
-      hsize_t(tmp.extent(0)), 
-      hsize_t(tmp.extent(1)), 
-      hsize_t(srfc ? tmp_srfc.extent(2) : tmp.extent(2)) 
+      hsize_t(arr.extent(0)), 
+      hsize_t(arr.extent(1)), 
+      hsize_t(srfc ? tmp_srfc.extent(2) : arr.extent(2)) 
     };
-    this->h5d.read(srfc ? tmp_srfc.data() : tmp.data(), H5::PredType::NATIVE_FLOAT, H5::DataSpace(3, ext), h5s);
+    this->h5d.read(srfc ? tmp_srfc.data() : arr.data(), H5::PredType::NATIVE_FLOAT, H5::DataSpace(3, ext), h5s);
 
-    return blitz::safeToReturn((srfc ? tmp_srfc : tmp) + 0);
+    return blitz::safeToReturn((srfc ? tmp_srfc : arr) + 0);
   }
 
   auto h5load_timestep(
@@ -235,6 +237,44 @@ class Plotter_t<3> : public PlotterCommon
     // other dataset are of the size x*z, resize tmp
     tmp.resize(n[0]-1, n[1]-1, n[2]-1);
     tmp_srfc.resize(n[0]-1, n[1]-1, 1);
+
+    // init refined data
+    try
+    {
+      this->h5f.openDataSet("X refined").getSpace().getSimpleExtentDims(n, NULL); 
+      this->map["refined x"] = n[0]-1;
+      this->map["refined y"] = n[1]-1;
+      this->map["refined z"] = n[2]-1;
+      tmp_ref.resize(n[0], n[1], n[2]);
+      h5load(file + "/const.h5", "X refined");
+      this->map["refined dx"] = tmp_ref(1,0,0) - tmp_ref(0,0,0);
+      h5load(file + "/const.h5", "Y refined");
+      this->map["refined dy"] = tmp_ref(0,1,0) - tmp_ref(0,0,0);
+      h5load(file + "/const.h5", "Z refined");
+      this->map["refined dz"] = tmp_ref(0,0,1) - tmp_ref(0,0,0);
+      this->CellVol_ref = this->map["refined dx"] * this->map["refined dy"] * this->map["refined dz"];
+      tmp_ref.resize(n[0]-1, n[1]-1, n[2]-1);
+    }
+    catch(...) // for pre-refinement simulations that didnt store refined stuff
+    {
+      this->map["refined x"] = this->map["x"]; 
+      this->map["refined y"] = this->map["y"]; 
+      this->map["refined z"] = this->map["z"]; 
+      this->map["refined dx"] = this->map["dx"]; 
+      this->map["refined dy"] = this->map["dy"]; 
+      this->map["refined dz"] = this->map["dz"]; 
+      this->CellVol_ref = this->CellVol;
+      tmp_ref.resize(tmp.shape());
+    }
+
+    for (auto const& x : this->map)
+{
+    std::cout << x.first  // string (key)
+              << ':'
+              << x.second // string's value
+              << std::endl;
+}
+
   }
 };
 
diff --git a/UWLCM_plotters/include/PlotterCommon.hpp b/UWLCM_plotters/include/PlotterCommon.hpp
index 1c75c83..5501425 100644
--- a/UWLCM_plotters/include/PlotterCommon.hpp
+++ b/UWLCM_plotters/include/PlotterCommon.hpp
@@ -14,7 +14,7 @@ class PlotterCommon
   std::map<std::string, double> map;
   std::map<std::string, arr_prof_t> map_prof;
   blitz::Array<float, 1> timesteps;
-  double CellVol, DomainSurf, DomainVol;
+  double CellVol, DomainSurf, DomainVol, CellVol_ref;
 
   protected:
   H5::H5File h5f;
@@ -61,6 +61,7 @@ class PlotterCommon
     auto attr = h5g.openAttribute(attr_name);
     notice_macro(std::string("about to read attribute value"))
     attr.read(attr.getDataType(), &ret);
+    notice_macro(std::string("attribute value read: ") + std::to_string(ret))
     return ret;
   }
 
@@ -120,6 +121,20 @@ class PlotterCommon
       h5s.getSimpleExtentDims(&n, NULL);
       map_prof.emplace("rhod", arr_prof_t(n));
       h5d.read(map_prof["rhod"].data(), H5::PredType::NATIVE_FLOAT);
+
+      // read dry air density profile on refined grid
+      try
+      {
+        h5load(file + "/const.h5", "refined rhod");
+        h5s.getSimpleExtentDims(&n, NULL);
+        map_prof.emplace("refined rhod", arr_prof_t(n));
+        h5d.read(map_prof["refined rhod"].data(), H5::PredType::NATIVE_FLOAT);
+      }
+      catch(...) // for pre-refinement simulations, use rhod as refined rhod
+      {
+        map_prof.emplace("refined rhod", map_prof["rhod"].copy());
+	std::cerr << "refined rhod as copy of rhod: " << map_prof["refined rhod"];
+      }
     }
   }
 };
diff --git a/UWLCM_plotters/include/PlotterMicro.hpp b/UWLCM_plotters/include/PlotterMicro.hpp
index 68dd71a..e4a2aa4 100644
--- a/UWLCM_plotters/include/PlotterMicro.hpp
+++ b/UWLCM_plotters/include/PlotterMicro.hpp
@@ -20,57 +20,81 @@ class PlotterMicro_t : public Plotter_t<NDims>
   const double L_evap = 2264.76e3; // latent heat of evaporation [J/kg]
 
   public:
+  void multiply_by_rhod(arr_t &arr)
+  {
+    if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0))
+      arr *= this->map_prof["rhod"](this->LastIndex);
+    else if(arr.extent(NDims-1) == this->map_prof["refined rhod"].extent(0))
+      arr *= this->map_prof["refined rhod"](this->LastIndex);
+    else
+      throw std::runtime_error("multiply_by_rhod: input array is neither normal grid size nor refined grid size");
+  }
+
+  void multiply_by_CellVol(arr_t &arr)
+  {
+    if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0))
+      arr *= this->CellVol;
+    else if(arr.extent(NDims-1) == this->map_prof["refined rhod"].extent(0))
+      arr *= this->CellVol_ref;
+    else
+      throw std::runtime_error("multiply_by_CellVol: input array is neither normal grid size nor refined grid size");
+  }
+
   // functions for diagnosing fields
   //
   // aerosol droplets mixing ratio
   auto h5load_ra_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("aerosol_rw_mom3", at) * 4./3. * 3.1416 * 1e3;
+      return arr_t(this->h5load_timestep("aerosol_rw_mom3", at) * 4./3. * 3.1416 * 1e3);
+    
     else if(this->micro == "blk_1m")
+    {
       res = 0;
-    return blitz::safeToReturn(res + 0);
+      return res;
+     // return blitz::safeToReturn(res + 0);
+    }
   }
 
   // cloud droplets mixing ratio
   auto h5load_rc_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3;
+      return arr_t(this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3);
     else if(this->micro == "blk_1m")
-      res = this->h5load_timestep("rc", at);
+      return arr_t(this->h5load_timestep("rc", at));
     else if(this->micro == "blk_2m")
-      res = this->h5load_timestep("rc", at);
-    return blitz::safeToReturn(res + 0);
+      return arr_t(this->h5load_timestep("rc", at));
   }
 
   // rain droplets mixing ratio
   auto h5load_rr_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3;
+      return arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3);
     else if(this->micro == "blk_1m")
-      res = this->h5load_timestep("rr", at);
+      return arr_t(this->h5load_timestep("rr", at));
     else if(this->micro == "blk_2m")
-      res = this->h5load_timestep("rr", at);
-    return blitz::safeToReturn(res + 0);
+      return arr_t(this->h5load_timestep("rr", at));
   }
 
   // activated drops mixing ratio
   auto h5load_ract_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) 
   {
     if(this->micro == "lgrngn")
     {
-      res = this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3;
-      res += arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3);
+      return arr_t(
+	       arr_t(this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3) + 
+               arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3)
+	     );
     }
     else if(this->micro == "blk_1m")
     {
@@ -82,34 +106,38 @@ class PlotterMicro_t : public Plotter_t<NDims>
       res = this->h5load_timestep("rc", at);
       res += arr_t(this->h5load_timestep("rr", at));
     }
-    return blitz::safeToReturn(res + 0);
+   // return blitz::safeToReturn(res + 0);
+    return res;
   }
 
   // cloud droplets concentration [1/kg]
   auto h5load_nc_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("cloud_rw_mom0", at);
+      return arr_t(this->h5load_timestep("cloud_rw_mom0", at));
     else if(this->micro == "blk_1m")
+    {
       res = 0;
+      return res;
+     // return blitz::safeToReturn(res + 0);
+    }
     else if(this->micro == "blk_2m")
-      res = this->h5load_timestep("nc", at);
-    return blitz::safeToReturn(res + 0);
+      return arr_t(this->h5load_timestep("nc", at));
   }
 
   // precipitation flux [W/m2]
   auto h5load_prflux_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  )// -> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
     {
-      res = this->h5load_timestep("precip_rate", at)
+      return arr_t(this->h5load_timestep("precip_rate", at)
               *  4./3 * 3.14 * 1e3 // to get mass
               / this->CellVol    // averaged over cell volume, TODO: make precip rate return specific moment? wouldnt need the dx and dy
-              * L_evap;
+              * L_evap);
     }
     else if(this->micro == "blk_1m")
       try
@@ -125,19 +153,21 @@ class PlotterMicro_t : public Plotter_t<NDims>
       {
         res = 0;
       }
-    return blitz::safeToReturn(res + 0);
+   // return blitz::safeToReturn(res + 0);
+     return res;
   }
 
   // RH
   auto h5load_RH_timestep(
     int at
-  ) -> decltype(blitz::safeToReturn(arr_t() + 0))
+  ) //-> decltype(blitz::safeToReturn(arr_t() + 0))
   {
     if(this->micro == "lgrngn")
-      res = this->h5load_timestep("RH", at);
+      return arr_t(this->h5load_timestep("RH", at));
     else if(this->micro == "blk_1m")
       res = 0;
-    return blitz::safeToReturn(res + 0);
+   // return blitz::safeToReturn(res + 0);
+     return res;
   }
 
 
diff --git a/drawbicyc/include/cases/RICO11/plots.hpp b/drawbicyc/include/cases/RICO11/plots.hpp
index e841b4f..0a4ccd3 100644
--- a/drawbicyc/include/cases/RICO11/plots.hpp
+++ b/drawbicyc/include/cases/RICO11/plots.hpp
@@ -1,14 +1,6 @@
 #pragma once
 
 const std::vector<std::string> series_rico({
-// "clb_bigrain_mean_rd",
-// "clb_bigrain_mean_kappa",
-// "clb_bigrain_mean_conc",
-// "clb_bigrain_mean_inclt",
-//, "clb_bigrain_mean_gccn_fraction"
-
-//"cl_acnv25_rico",
-//"cl_accr25_rico",
 "RH_max",
 "cloud_cover_rico",
 "min_cloud_base_rico",
@@ -17,14 +9,25 @@ const std::vector<std::string> series_rico({
 "lwp",
 "rwp",
 "surf_precip",
-"acc_precip"
-,"nc"
-,"nr"
-,"cl_nc_rico"
-,"cl_nr_rico"
-,"cloud_avg_supersat"
-,"wvarmax"
-,"cl_meanr" //TODO: zmienic maske na rico
+"acc_precip",
+"cl_nc_rico",
+"cl_nr_rico",
+"cloud_avg_supersat",
+"wvarmax",
+"cl_meanr", //TODO: zmienic maske na rico
+"sd_conc"
+
+
+// "clb_bigrain_mean_rd",
+// "clb_bigrain_mean_kappa",
+// "clb_bigrain_mean_conc",
+// "clb_bigrain_mean_inclt",
+//, "clb_bigrain_mean_gccn_fraction"
+
+//"cl_acnv25_rico",
+//"cl_accr25_rico",
+//,"nc"
+//,"nr"
 //TODO (po usprawnieniu cloud mask i ujednoliceniu tego:
 /*
 ,"cl_avg_cloud_rad"
@@ -34,7 +37,6 @@ const std::vector<std::string> series_rico({
 //,"rd_geq_0.8um_conc"
 //,"cl_rd_lt_0.8um_conc"
 //,"cl_rd_geq_0.8um_conc"
-,"sd_conc"
 
 /*
  "cloud_base",
diff --git a/drawbicyc/include/gnuplot_series_set_labels.hpp b/drawbicyc/include/gnuplot_series_set_labels.hpp
index 49757f1..edc4c02 100644
--- a/drawbicyc/include/gnuplot_series_set_labels.hpp
+++ b/drawbicyc/include/gnuplot_series_set_labels.hpp
@@ -151,7 +151,7 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt)
     gp << "set title 'average cloud drop conc [1/cm^3]'\n";
   else if (plt == "ntot")
     gp << "set title 'average particle conc [1/cm^3]'\n";
-  else if (plt == "cl_nc")
+  else if (plt == "cl_nc" || plt == "cl_nc_rico")
   {
     gp << "set title 'average cloud drop conc [1/cm^3] in cloudy cells'\n";
     gp << "set xlabel ''\n";
@@ -163,7 +163,7 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt)
     gp << "set xlabel ''\n";
     gp << "set ylabel ''\n";
   }
-  else if (plt == "cl_nr")
+  else if (plt == "cl_nr" || plt == "cl_nr_rico")
   {
     gp << "set title 'average rain drop conc [1/cm^3] in cloudy cells'\n";
     gp << "set xlabel ''\n";
diff --git a/drawbicyc/include/plot_series.hpp b/drawbicyc/include/plot_series.hpp
index 2167e77..c4dffa2 100644
--- a/drawbicyc/include/plot_series.hpp
+++ b/drawbicyc/include/plot_series.hpp
@@ -129,8 +129,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg
           typename Plotter_t::arr_t snap(tmp); 
           snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg
-          snap *= rhod; // water per cubic metre (should be wet density...)
-          plotter.k_i = blitz::sum(snap, plotter.LastIndex) * n["dz"]; // LWP [g/m2] in the column 
+          plotter.multiply_by_rhod(snap); 
+          plotter.k_i = blitz::sum(snap, plotter.LastIndex) * n["refined dz"]; // LWP [g/m2] in the column 
           plotter.k_i = where(plotter.k_i > 20 , 1 , 0); // cloudiness as in Ackermann et al. 
           res_series[plt](at) = blitz::mean(plotter.k_i);
         }
@@ -245,19 +245,19 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           {
             auto tmp = plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"]) * 4./3. * 3.1416 * 1e3;
             typename Plotter_t::arr_t snap(tmp);
-            snap *= rhod;
+            plotter.multiply_by_rhod(snap); 
             res_series[plt](at) += blitz::mean(snap);
           }
           {
             auto tmp = plotter.h5load_timestep("rain_rw_mom3", at * n["outfreq"]) * 4./3. * 3.1416 * 1e3;
             typename Plotter_t::arr_t snap(tmp);
-            snap *= rhod;
+            plotter.multiply_by_rhod(snap); 
             res_series[plt](at) += blitz::mean(snap);
           }
           {
             auto tmp = plotter.h5load_timestep("rv", at * n["outfreq"]);
             typename Plotter_t::arr_t snap(tmp);
-            snap *= rhod;
+            plotter.multiply_by_rhod(snap); 
             res_series[plt](at) += blitz::mean(snap);
           } 
         }
@@ -342,7 +342,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
             auto tmp2 = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]);
             typename Plotter_t::arr_t snap_mom(tmp2);
             com_N_c(at) = snap_mom(com_x_idx(at), com_z_idx(at)); // 0th raw moment / mass [1/kg]
-            snap_mom *= rhod; // now per m^3
+            plotter.multiply_by_rhod(snap); 
             res_series[plt](at) = snap_mom(com_x_idx(at), com_z_idx(at));
           }
           else 
@@ -457,7 +457,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap(tmp);
           snap /= 1e6; // per cm^3
-          snap *= rhod; // b4 it was per milligram
+          plotter.multiply_by_rhod(snap); 
           res_series[plt](at) = blitz::mean(snap); 
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
@@ -470,7 +470,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           auto tmp = plotter.h5load_timestep("all_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap(tmp);
           snap /= 1e6; // per cm^3
-          snap *= rhod; // b4 it was per milligram
+          plotter.multiply_by_rhod(snap); 
           res_series[plt](at) = blitz::mean(snap); 
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
@@ -483,7 +483,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           auto tmp = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap(tmp);
           snap /= 1e6; // per cm^3
-          snap *= rhod; // b4 it was per milligram
+          plotter.multiply_by_rhod(snap); 
           res_series[plt](at) = blitz::mean(snap); 
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
@@ -496,7 +496,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           // cloud fraction (cloudy if N_c > 20/cm^3)
           auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap(tmp);
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           typename Plotter_t::arr_t snap2;
           snap2.resize(snap.shape());
@@ -522,7 +522,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux
 
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap2 *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap2); 
           snap2 /= 1e6; // per cm^3
 
           snap2 *= snap;
@@ -541,12 +541,12 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           // cloud fraction (cloudy if N_c > 20/cm^3)
           auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap(tmp);
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           auto tmp2 = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap2(tmp2);
-          snap2 *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap2); 
           snap2 /= 1e6; // per cm^3
           snap2 *= snap;
           if(blitz::sum(snap) > 0)
@@ -569,7 +569,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
 
           auto tmp2 = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap2(tmp2);
-          snap2 *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap2); 
           snap2 /= 1e6; // per cm^3
 
           snap2 *= snap;
@@ -599,7 +599,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           // cloud fraction (cloudy if N_c > 20/cm^3)
           auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap(tmp);
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux
@@ -609,7 +609,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           cloudy_column = where(cloudy_column > 0, 1, 0);
           plotter.k_i = where(cloudy_column == 0, 0, plotter.k_i);
           if(blitz::sum(cloudy_column) > 0)
-            res_series[plt](at) = double(blitz::sum(plotter.k_i)) / blitz::sum(cloudy_column) * n["dz"];
+            res_series[plt](at) = double(blitz::sum(plotter.k_i)) / blitz::sum(cloudy_column) * n["refined dz"];
           else
             res_series[plt](at) = 0;
         }
@@ -624,7 +624,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           // cloud fraction (cloudy if N_c > 20/cm^3)
           auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap(tmp);
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux
@@ -665,7 +665,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           cloudy_column = where(cloudy_column > 0, 1, 0);
           plotter.k_i = where(cloudy_column == 0, 1e6, plotter.k_i); // 1e6 denotes no clouds in the column
           if(blitz::sum(cloudy_column) > 0)
-            res_series[plt](at) = blitz::min(plotter.k_i) * n["dz"];
+            res_series[plt](at) = blitz::min(plotter.k_i) * n["refined dz"];
           else
             res_series[plt](at) = 0;
         }
@@ -723,7 +723,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           auto tmp = plotter.h5load_timestep("rd_rng000_mom3", at * n["outfreq"]) * 4./3. * 3.14 * rho_dry * 1e3;
           typename Plotter_t::arr_t snap(tmp);
-          snap *= rhod * plotter.CellVol; // turn mixing ratio in g/kg to total mass in g
+          plotter.multiply_by_rhod(snap); 
+          plotter.multiply_by_CellVol(snap); 
           res_series[plt](at) = blitz::sum(snap); 
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
@@ -765,7 +766,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
           typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
 
@@ -791,7 +792,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
           typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
 
@@ -861,11 +862,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         try
         {
           {
-            auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg
-            typename Plotter_t::arr_t snap(tmp); 
-            snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg
-            snap *= rhod; // water per cubic metre (should be wet density...)
-            res_series[plt](at) = blitz::mean(snap); 
+            typename Plotter_t::arr_t rl(plotter.h5load_rc_timestep(at * n["outfreq"]));
+            typename Plotter_t::arr_t rr(plotter.h5load_rr_timestep(at * n["outfreq"]));
+	    rl = (rl + rr) * 1e3; // g/kg
+            plotter.multiply_by_rhod(rl); 
+            res_series[plt](at) = blitz::mean(rl); 
           }
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
@@ -877,7 +878,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           {
             typename Plotter_t::arr_t snap(plotter.h5load_rr_timestep(at * n["outfreq"]));
-            snap *= rhod * 1e3; // water per cubic metre (should be wet density...) & g/kg
+            snap *= 1e3;
+            plotter.multiply_by_rhod(snap); 
             res_series[plt](at) = blitz::mean(snap); 
           }
         }
@@ -1172,11 +1174,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
           typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"]));
-          snap2 *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap2); 
           snap2 /= 1e6; // per cm^3
           snap2 *= snap;
           if(blitz::sum(snap) > 0)
@@ -1193,11 +1195,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
           typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"]));
-          snap2 *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap2); 
           snap2 /= 1e6; // per cm^3
           snap2 *= snap;
           if(blitz::sum(snap) > 0)
@@ -1213,7 +1215,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
           typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"]));
@@ -1236,7 +1238,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
           typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           typename Plotter_t::arr_t snap_m0(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"]));
@@ -1258,7 +1260,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"]));
           snap2 /= 1e6; // per cm^3
-          snap2 *= rhod; // b4 it was per milligram
+          plotter.multiply_by_rhod(snap2); 
           res_series[plt](at) = blitz::mean(snap2); 
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
@@ -1270,7 +1272,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"]));
           snap2 /= 1e6; // per cm^3
-          snap2 *= rhod; // b4 it was per milligram
+          plotter.multiply_by_rhod(snap2); 
           res_series[plt](at) = blitz::mean(snap2); 
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
@@ -1284,11 +1286,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
           typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_geq_0.8um_rw_mom0", at * n["outfreq"]));
-          snap2 *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap2); 
           snap2 /= 1e6; // per cm^3
           snap2 *= snap;
           if(blitz::sum(snap) > 0)
@@ -1305,7 +1307,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_geq_0.8um_rw_mom0", at * n["outfreq"]));
           snap2 /= 1e6; // per cm^3
-          snap2 *= rhod; // b4 it was per milligram
+          plotter.multiply_by_rhod(snap2); 
           res_series[plt](at) = blitz::mean(snap2); 
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
@@ -1318,11 +1320,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
           typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_lt_0.8um_rw_mom0", at * n["outfreq"]));
-          snap2 *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap2); 
           snap2 /= 1e6; // per cm^3
           snap2 *= snap;
           if(blitz::sum(snap) > 0)
@@ -1338,7 +1340,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_lt_0.8um_rw_mom0", at * n["outfreq"]));
           snap2 /= 1e6; // per cm^3
-          snap2 *= rhod; // b4 it was per milligram
+          plotter.multiply_by_rhod(snap2); 
           res_series[plt](at) = blitz::mean(snap2); 
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
@@ -1351,7 +1353,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
           typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
-          snap *= rhod; // b4 it was specific moment
+          plotter.multiply_by_rhod(snap); 
           snap /= 1e6; // per cm^3
           snap = iscloudy(snap); // cloudiness mask
           typename Plotter_t::arr_t snap_m0(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
@@ -1933,11 +1935,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
     }
     else if (plt == "lwp")
     {
-      res_series[plt] *= (n["z"] - 1) * n["dz"]; // top and bottom cells are smaller
+      res_series[plt] *= (n["refined z"] - 1) * n["refined dz"]; // top and bottom cells are smaller
     }
     else if (plt == "rwp")
     {
-      res_series[plt] *= (n["z"] - 1) * n["dz"]; // top and bottom cells are smaller
+      res_series[plt] *= (n["refined z"] - 1) * n["refined dz"]; // top and bottom cells are smaller
     }
     else if (plt == "cwp")
     {