diff --git a/UWLCM_plotters/include/Plotter2d.hpp b/UWLCM_plotters/include/Plotter2d.hpp
index 2bed638..63ea185 100644
--- a/UWLCM_plotters/include/Plotter2d.hpp
+++ b/UWLCM_plotters/include/Plotter2d.hpp
@@ -198,6 +198,17 @@ class Plotter_t<2> : public PlotterCommon
     // edge cells are smaller
     dv(blitz::Range(0,0),blitz::Range::all()) /= 2.;
     dv(blitz::Range::all(),blitz::Range(0,0)) /= 2.;
+    dv(blitz::Range::all(),blitz::Range(dv.cols() - 1, dv.cols() - 1)) /= 2.;
+    dv(blitz::Range(dv.rows() - 1, dv.rows() - 1), blitz::Range::all()) /= 2.;
+    
+    // Printing the modified dv array
+    //std::cout << "Modified dv array:\n";
+    //for (int i = 0; i < dv.extent(0); ++i) {
+    //   for (int j = 0; j < dv.extent(1); ++j) {
+    //        std::cout << dv(i, j) << " ";
+    //   }
+    //   std::cout << "\n";
+    //}
 
     // other dataset are of the size x*z, resize tmp
     tmp.resize(n[0]-1, n[1]-1);
diff --git a/UWLCM_plotters/include/PlotterCommon.hpp b/UWLCM_plotters/include/PlotterCommon.hpp
index 1c75c83..9d33f4b 100644
--- a/UWLCM_plotters/include/PlotterCommon.hpp
+++ b/UWLCM_plotters/include/PlotterCommon.hpp
@@ -30,10 +30,10 @@ class PlotterCommon
   {
     if(h5f.getFileName() != file)
     {
-      notice_macro("about to close current file: " << h5f.getFileName())
+      //notice_macro("Teraz about to close current file: " << h5f.getFileName())
       h5f.close();
 
-      notice_macro("about to open file: " << file)
+      //notice_macro("Teraaz 1about to open file: " << file)
       h5f.openFile(file, H5F_ACC_RDONLY);
     }
 
@@ -46,10 +46,10 @@ class PlotterCommon
   {
     if(h5f.getFileName() != file)
     {
-      notice_macro("about to close current file: " << h5f.getFileName())
+      //notice_macro("Teraz 2 about to close current file: " << h5f.getFileName())
       h5f.close();
   
-      notice_macro("about to open file: " << file)
+      //notice_macro("Teraz 2about to open file: " << file)
       h5f.openFile(file, H5F_ACC_RDONLY);
     }
   
@@ -85,7 +85,7 @@ class PlotterCommon
     file(file)
   {
     // init h5f
-    notice_macro("about to open file: " << file << "/const.h5")
+    //notice_macro("Teraz 3about to open file: " << file << "/const.h5")
     h5f.openFile(file + "/const.h5", H5F_ACC_RDONLY);
 
     // init dt and outfreq
diff --git a/UWLCM_plotters/include/PlotterMicro.hpp b/UWLCM_plotters/include/PlotterMicro.hpp
index 1945a10..1da8f28 100644
--- a/UWLCM_plotters/include/PlotterMicro.hpp
+++ b/UWLCM_plotters/include/PlotterMicro.hpp
@@ -106,10 +106,16 @@ class PlotterMicro_t : public Plotter_t<NDims>
   {
     if(this->micro == "lgrngn")
     {
+
+      //std::cout << "CellVol = " << this->CellVol << std::endl;
+      //std::cout << "dv = " << this->dv << std::endl;
+      //std::cerr << this->dv;
       res = 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;
+            *  4./3 * 3.14 * 1e3 // to get mass
+              //5000.
+	    / this->dv    // averaged over cell volume, TODO: make precip rate return specific moment? wouldnt need the dx and dy
+             // this->CellVol    // averaged over cell volume, TODO: make precip rate return specific moment? wouldnt need the dx and dy
+            * L_evap;
     }
     else if(this->micro == "blk_1m")
       try
@@ -388,7 +394,14 @@ class PlotterMicro_t : public Plotter_t<NDims>
                      / 1e3         // to m^3 of water
                      * 1e3;        // to mm
   }
-
+  // accumulated volume precipitation [m^3]
+  double calc_acc_surf_precip_volume(double prec_vol)
+  {
+    if(this->micro == "lgrngn")
+      return calc_acc_surf_precip(prec_vol) * this->DomainSurf / 1000.;
+    if(this->micro == "blk_1m")
+      return calc_acc_surf_precip(prec_vol) * this->DomainSurf / 1000.;// to m^3 of water
+  }
   // droplet removal rate (at boundaries) since last output [1/(cm^3 s)]
   double calc_prtcl_removal(double prtcl_removal_diff)
   {
diff --git a/cases/Lasher_Trapp/Lasher_Trapp_series_comparison.py b/cases/Cumulus_Congestus/Cumulus_Congestus_series_comparison.py
similarity index 86%
rename from cases/Lasher_Trapp/Lasher_Trapp_series_comparison.py
rename to cases/Cumulus_Congestus/Cumulus_Congestus_series_comparison.py
index 876bf10..40ac10b 100644
--- a/cases/Lasher_Trapp/Lasher_Trapp_series_comparison.py
+++ b/cases/Cumulus_Congestus/Cumulus_Congestus_series_comparison.py
@@ -12,23 +12,23 @@
 # activate latex text rendering
 rc('text', usetex=True)
 
-Lasher_Trapp_vars = ["lwp", "rwp", "surf_precip", "acc_precip", "cl_nc"]
+Cumulus_Congestus_vars = ["RH_max", "cloud_top_height", "surf_precip", "acc_precip", "acc_vol_precip"]
 
 # init the plot
 nplotx = 2
 nploty= 3
 fig, axarr = plt.subplots(nplotx,nploty)
 
-plot_series( Lasher_Trapp_vars, 0, nplotx, nploty, axarr, xscaledict, yscaledict, xlimdict_series, ylimdict_series, xlabel='Time [h]')
+plot_series( Cumulus_Congestus_vars, 0, nplotx, nploty, axarr, xscaledict, yscaledict, xlimdict_series, ylimdict_series, xlabel='Time [h]')
 
 # legend font size
 plt.rcParams.update({'font.size': 8})
 
 # hide axes on empty plots
-if len( Lasher_Trapp_vars) % nploty == 0:
+if len( Cumulus_Congestus_vars) % nploty == 0:
   nemptyplots = 0
 else:
-  nemptyplots = nploty - len( Lasher_Trapp_vars) % nploty
+  nemptyplots = nploty - len( Cumulus_Congestus_vars) % nploty
   emptyplots = np.arange(nploty - nemptyplots, nploty)
 for empty in emptyplots:
   axarr[nplotx-1, empty].axis('off')
diff --git a/cases/Lasher_Trapp/plot_ranges.py b/cases/Cumulus_Congestus/plot_ranges.py
similarity index 77%
rename from cases/Lasher_Trapp/plot_ranges.py
rename to cases/Cumulus_Congestus/plot_ranges.py
index 6dc3ccc..f4e2d51 100644
--- a/cases/Lasher_Trapp/plot_ranges.py
+++ b/cases/Cumulus_Congestus/plot_ranges.py
@@ -10,6 +10,10 @@
   "sat_RH" : "linear",
   "rad_flx" : "linear",
   "lwp" : "linear",
+  "RH_max" : "linear",
+  "cloud_top_height" : "linear",
+  "total_droplets_number" : "linear",
+  "acc_precip_vol" : "linear",
   "rwp" : "linear",
   "er" : "linear",
   "wvarmax" : "linear",
@@ -34,6 +38,10 @@
   "sat_RH" : "linear",
   "rad_flx" : "linear",
   "lwp" : "linear",
+  "RH_max" : "linear",
+  "cloud_top_height" : "linear",
+  "total_droplets_number" : "linear",
+  "acc_precip_vol" : "linear",
   "rwp" : "linear",
   "er" : "linear",
   "wvarmax" : "linear",
@@ -56,6 +64,8 @@
   "wvar" : None,
   "w3rd" : None,
   "sat_RH" : None,
+  #"RH_max" : None,
+  #"cloud_top_height" : None,
   "rad_flx" : None,
   "gccn_rw_cl" : (0,90),
   "non_gccn_rw_cl" : (0,12),
@@ -71,6 +81,8 @@
   "wvar" : (0,3000),
   "w3rd" : (0,3000),
   "sat_RH" : (0,3000),
+  "RH_max" : (0,3000),
+  "cloud_top_height" : (0,3000),
   "rad_flx" : (0,3000),
   "gccn_rw_cl" : (0,3000),
   "non_gccn_rw_cl" : (0,3000),
@@ -81,24 +93,32 @@
   "clfrac" :  None,
   "cl_nc" :  None,
   "lwp" :  None,
+  "RH_max" : None,
+  "cloud_top_height" : None,
+  "total_droplets_number" : None,
   "rwp" :  None,
   "er" :  None,
   "wvarmax" :  None,
   "surf_precip" :  None,
   "acc_precip" :  None,
+  "acc_precip_vol" : None,
   "cl_gccn_conc" :  None,
-  "cloud_base" :  None
+  "cloud_base" :  None  
   }
 
 ylimdict_series = {
   "clfrac" : None,
   "cl_nc" : None,
   "lwp" : None,
+  "RH_max" : None,
+  "cloud_top_height" : None,
+  "total_droplets_number" : None,
+  "acc_precip_vol" : None,
   "rwp" : None,
   "er" : None,
   "wvarmax" : None,
   "surf_precip" : None,
   "acc_precip" : None,#(0,0.07),
   "cl_gccn_conc" : (1e-6, 1),
-  "cloud_base" : None
-  }
+  "cloud_base" : None  
+}
diff --git a/drawbicyc/average.cpp b/drawbicyc/average.cpp
index c5cd768..5dca5f2 100644
--- a/drawbicyc/average.cpp
+++ b/drawbicyc/average.cpp
@@ -193,7 +193,7 @@ int main(int argc, char* argv[])
  const string profiles_suffix = string("_profiles_")+argv[2]+string("_")+argv[3];
  string plot_type;
  // determine type of plots based on the name of the first file
- const string types[] = {"rico", "dycoms", "moist_thermal"};
+ const string types[] = {"rico", "dycoms", "moist_thermal", "cumulus_congestus"};
  for(auto type : types)
  {
    if(hasEnding(string(argv[4]), type))
diff --git a/drawbicyc/drawbicyc.cpp b/drawbicyc/drawbicyc.cpp
index fd01288..a12306f 100644
--- a/drawbicyc/drawbicyc.cpp
+++ b/drawbicyc/drawbicyc.cpp
@@ -19,7 +19,7 @@ int main(int argc, char** argv)
     ("qv_qc_2_6_10_min", po::value<bool>()->default_value(false) , "plot comparison of qv and qc fields at 2, 6 and 10 min?")
     ("dir", po::value<std::string>()->required() , "directory containing out_lgrngn")
     ("micro", po::value<std::string>()->required(), "one of: blk_1m, blk_2m, lgrngn")
-    ("type", po::value<std::string>()->required(), "one of: dycoms, moist_thermal, rico, Lasher_Trapp, gccn_ccn_conc")//, base_prflux_vs_clhght")
+    ("type", po::value<std::string>()->required(), "one of: dycoms, moist_thermal, rico, cumulus_congestus,  gccn_ccn_conc")//, base_prflux_vs_clhght")
   ;
 
   po::variables_map vm;
@@ -33,8 +33,8 @@ int main(int argc, char** argv)
 
   // handling the "type" option
   std::string type = vm["type"].as<std::string>();
-  if(type != "dycoms" && type != "moist_thermal" && type != "rico" && type != "pi_chamber" && type != "Lasher_Trapp" && type != "pi_chamber_icmw" && type != "gccn_ccn_conc")
-    throw std::runtime_error("Unrecognized 'type' option, only dycoms, rico, moist_thermal, pi_chamber, pi_chamber_icmw, Lasher_Trapp, gccn_ccn_conc  available now");//, base_prflux_vs_clhght available now");
+  if(type != "dycoms" && type != "moist_thermal" && type != "rico" && type != "pi_chamber" && type !=  "cumulus_congestus" && type != "pi_chamber_icmw" && type != "gccn_ccn_conc")
+    throw std::runtime_error("Unrecognized 'type' option, only dycoms, rico, moist_thermal, pi_chamber, pi_chamber_icmw, cumulus_congestus, gccn_ccn_conc  available now");//, base_prflux_vs_clhght available now");
 
   // should profiles be normalized by inversion height
   const bool normalize_prof = type == "dycoms";
@@ -56,9 +56,10 @@ int main(int argc, char** argv)
   H5::DataSet h5d = h5f.openDataSet("G");
   H5::DataSpace h5s = h5d.getSpace();
   int NDims = h5s.getSimpleExtentNdims();
-  
-  // detecting if subgrid model was on
+  std::cout << NDims <<std::endl;
+  //detecting if subgrid model was on
   bool sgs = true;
+  std::cout<<"WARTOSC SGS"<<std::endl;
   try 
   {
     auto h5g = h5f.openGroup("sgs");
@@ -69,26 +70,27 @@ int main(int argc, char** argv)
   }
 
   Plots plots(type, sgs);
+  //Plots plots(type, false);
 
   if(NDims == 2)
   {
 
     if(flag_series)   plot_series(PlotterMicro_t<2>(h5, micro), plots, type);
     if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof);
-//    if(flag_fields)   plot_fields(PlotterMicro_t<2>(h5, micro), plots, type);
-//    if(flag_qv_qc_2_6_10_min)   plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro));
+    //if(flag_fields)   plot_fields(PlotterMicro_t<2>(h5, micro), plots, type);
+    //if(flag_qv_qc_2_6_10_min)   plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro));
   }
-  else if(NDims == 3)
-  {
-    if(flag_series)   plot_series(PlotterMicro_t<3>(h5, micro), plots, type);
-    if(flag_profiles) plot_profiles(PlotterMicro_t<3>(h5, micro), plots, type, normalize_prof);
-//    if(flag_fields)   plot_fields(PlotterMicro_t<3>(h5, micro), plots, type);
-//    if(flag_qv_qc_2_6_10_min)   plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro));
+//  else if(NDims == 3)
+//  {
+//    if(flag_series)   plot_series(PlotterMicro_t<3>(h5, micro), plots, type);
+    //if(flag_profiles) plot_profiles(PlotterMicro_t<3>(h5, micro), plots, type, normalize_prof);
+    //if(flag_fields)   plot_fields(PlotterMicro_t<3>(h5, micro), plots, type);
+    //if(flag_qv_qc_2_6_10_min)   plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro));
 //    if(flag_lgrngn_spec) {
 //      plot_lgrngn_spec_positions(PlotterMicro_t<3>(h5, "lgrngn"));
 //      plot_lgrngn_spec(PlotterMicro_t<3>(h5, "lgrngn"));
 //    }
-  }
+//  }
   else
     assert(false && "need 2d or 3d input data");
 
diff --git a/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp
new file mode 100644
index 0000000..544e360
--- /dev/null
+++ b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp
@@ -0,0 +1,112 @@
+#pragma once
+
+
+const std::vector<std::string> series_Cumulus_Congestus({
+  //"clfrac",
+  "lwp",
+  "rwp",
+  "cwp",
+  "lwm",
+  "cwm",
+  "rwm",
+  "surf_precip",
+  "acc_precip",
+  "acc_vol_precip",
+  //"cl_nc",
+  //"cloud_base",
+  //"RH_max",
+  "cloud_top_height",
+  "total_droplets_number",
+  "cl_meanr",
+  "cl_nc",
+  "cl_nr",
+  "nc",
+  "nr",
+  "com_mom0",
+  "com_mom1",
+  "com_mom2",
+  "cloud_cover_rico", 
+  //"surf_flux_latent",
+  //"surf_flux_sensible",
+  //"sd_conc_avg",
+  //"mass_dry",
+  //"cl_gccn_conc",
+  //"gccn_conc",
+  //"cl_non_gccn_conc",
+  //"non_gccn_conc",
+  //"cl_gccn_to_non_gccn_conc_ratio"
+  //, "cl_gccn_meanr"
+  //,"cl_avg_cloud_rad"
+  // "sd_conc_std_dev",
+  // "tot_water"
+});
+ 
+std::vector<std::string> profs_Cumulus_Congestus({
+  "actrw_rw_cl_conc",
+  "rliq",
+  "actrw_rw_cl",
+  "actrw_mom1",
+  "actrw_mom2",
+  "disp_r",
+  "actrw_rw_cl_sigma",
+  "sigma_r",
+  "mean_r",
+  //"00rtot",
+  //"rliq",
+  //"thl",
+  //"wvar",
+  "prflux",
+  //"clfrac",
+  //"sd_conc",
+  //"cl_nc",
+  //"cl_nc_up",
+  //"w",
+  //"u",
+  //"v",
+  //"coal_tele",
+  //"base_prflux_vs_clhght",
+  //"non_gccn_rw_cl",
+  //"gccn_rw_cl,"
+  //, "N_c",
+ // "actrw_reff_cl",
+ // "ratio_mean_volume_r_to_eff_r_cubed",
+ // "cloud_water_std",
+  //"rain_water_std",
+ // "cloud_std_dev"
+  //,"vel_div"
+  //, "nc_up"
+  //,"sat_RH_up"
+  //, "act_conc_up"
+  //, "nc_down"
+});
+
+std::vector<std::string> fields_Cumulus_Congestus({
+  //"rl",
+  "mass_conc_rain_drops",
+  "mass_conc_cloud_droplets",
+  "num_conc_droplets",
+  "LWC",
+  //"mass_conc_cloud",
+  "ratio_mean_volume_r_to_eff_r_cubed",
+  "std_dev_droplet_size_dist",
+  //"rr",
+  "num_conc_rain_drops",
+  "effective_radius",
+  //"na",
+  //"th",
+  //"rv"
+  //"u",
+  //"w"
+  //"sd_conc",
+  //"r_dry",
+  //"RH",
+  //"supersat",
+  //"nc",
+  //"rr",
+  //"nr",
+  //"ef", "na",
+  //"th",
+  //"rv",
+  //"u",
+  //"w"
+});
diff --git a/drawbicyc/include/cases/Lasher_Trapp/plots.hpp b/drawbicyc/include/cases/Lasher_Trapp/plots.hpp
deleted file mode 100644
index 5ce831f..0000000
--- a/drawbicyc/include/cases/Lasher_Trapp/plots.hpp
+++ /dev/null
@@ -1,63 +0,0 @@
-#pragma once
-
-
-const std::vector<std::string> series_Lasher_Trapp({
-  "clfrac",
-  "lwp",
-  "rwp",
-  "surf_precip",
-  "acc_precip",
-  "cl_nc",
-  "cloud_base",
-  "surf_flux_latent",
-  "surf_flux_sensible",
-  "sd_conc_avg",
-  //"mass_dry",
-  "cl_gccn_conc",
-  "gccn_conc",
-  "cl_non_gccn_conc",
-  "non_gccn_conc",
-  "cl_gccn_to_non_gccn_conc_ratio"
-  //, "cl_gccn_meanr"
-  //,"cl_avg_cloud_rad"
-  // "sd_conc_std_dev",
-  // // "tot_water"
-});
- 
-std::vector<std::string> profs_Lasher_Trapp({
-  "00rtot",
-  "rliq",
-  "thl",
-  "wvar",
-  "prflux",
-  "clfrac",
-  "sd_conc",
-  "cl_nc",
-  "cl_nc_up",
-  "w",
-  "u",
-  "v",
-  "base_prflux_vs_clhght",
-  "non_gccn_rw_cl",
-  "gccn_rw_cl,"
-  //, "N_c",
-  //,"vel_div"
-  //, "nc_up"
-  //,"sat_RH_up"
-  //, "act_conc_up"
-  //, "nc_down"
-});
-
-std::vector<std::string> fields_Lasher_Trapp({
-  "rl",
-  "nc",
-  "rr",
-  "nr",
-  //"ef", "na",
-  "th",
-  "rv",
-  "u",
-  "w"
-  //"sd_conc",//, "r_dry",
-  //"RH", "supersat",
-});
diff --git a/drawbicyc/include/gnuplot_profs_set_labels.hpp b/drawbicyc/include/gnuplot_profs_set_labels.hpp
index f2cfd8e..eb83cdf 100644
--- a/drawbicyc/include/gnuplot_profs_set_labels.hpp
+++ b/drawbicyc/include/gnuplot_profs_set_labels.hpp
@@ -13,7 +13,7 @@ void gnuplot_profs_set_labels(Gnuplot &gp, std::string plt, const bool normalize
 
   if (plt == "rliq")
   {
-    gp << "set title 'liquid water [g/kg]'\n";
+    gp << "set title 'liquid water in cloudy cells [g/kg]'\n";
   }
   if (plt == "gccn_rw")
   {
@@ -162,6 +162,10 @@ void gnuplot_profs_set_labels(Gnuplot &gp, std::string plt, const bool normalize
   {
     gp << "set title 'precipitation flux [W/m^2]'\n";
   }
+  else if (plt == "coal_tele")
+  {
+    gp << "set title 'coal tele mass flux [W/m^2]'\n";
+  }
   else if (plt == "rad_flx")
   {
     gp << "set title 'radiative flux [W/m2]'\n";
@@ -214,4 +218,28 @@ void gnuplot_profs_set_labels(Gnuplot &gp, std::string plt, const bool normalize
   {
     gp << "set title 'mean cloud droplet radius [um]'\n";
   }
+    else if (plt == "actrw_mom1")
+  {
+    gp << "set title 'mean radius of actrw droplets [um]'\n";
+  }
+  else if (plt == "actrw_mom2")
+  {
+    gp << "set title 'variance of actrw droplet size distribution [um^2] * 1e3'\n";
+  }
+  else if (plt == "actrw_sd")
+  {
+    gp << "set title 'standard deviation of actrw droplet size distribution [um]'\n";
+  }
+  else if (plt == "actrw_rw_cl_conc")
+  {
+    gp << "set title 'concentration of actrw droplet in cloudy cells [cm^-3]'\n";
+  }
+  else if (plt == "actrw_rw_cl")
+  {
+    gp << "set title 'mean radius of actrw droplet in cloudy cells [um]'\n";
+  }
+    else if (plt == "disp_r")
+  {
+    gp << "set title 'relative dispersion of r_{drop} [1]'\n";
+  }
 }
diff --git a/drawbicyc/include/gnuplot_series_set_labels.hpp b/drawbicyc/include/gnuplot_series_set_labels.hpp
index b258c9d..49757f1 100644
--- a/drawbicyc/include/gnuplot_series_set_labels.hpp
+++ b/drawbicyc/include/gnuplot_series_set_labels.hpp
@@ -187,6 +187,24 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt)
     gp << "set xlabel ''\n";
     gp << "set ylabel ''\n";
   }
+  else if (plt == "cloud_top_height")
+  {
+    gp << "set title 'cloud top height [m]'\n";
+    gp << "set xlabel 'time [min]'\n";
+    gp << "set ylabel 'height [m]'\n";
+  }
+  else if (plt == "total_droplets_number")
+  {
+    gp << "set title 'Total droplets number [#]'\n";
+    gp << "set xlabel 'time [min]'\n";
+    gp << "set ylabel ''\n";
+  }
+  else if (plt == "acc_vol_precip")
+  {
+    gp << "set title 'accumulated volume precipitation [mm^3]'\n";
+    gp << "set xlabel 'time [h]'\n";
+    gp << "set ylabel ''\n";
+  }
   else if (plt == "mass_dry")
     gp << "set title 'total dry mass [g]'\n";
   else if (plt == "RH_max")
@@ -198,13 +216,35 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt)
   else if (plt == "lwp")
   {
     gp << "set title 'liquid water path [g / m^2]'\n";
-    gp << "set xlabel ''\n";
+    gp << "set xlabel 'time[h]'\n";
+    gp << "set ylabel ''\n";
+  }
+  else if (plt == "lwm")
+  {
+    gp << "set title 'liquid water mass [kg]'\n";
+    gp << "set xlabel 'time[h]'\n";
+    gp << "set ylabel ''\n";
+  }else if (plt == "cwm")
+  {
+    gp << "set title 'cloud water mass [kg]'\n";
+    gp << "set xlabel 'time[h]'\n";
+    gp << "set ylabel ''\n";
+  }else if (plt == "rwm")
+  {
+    gp << "set title 'rain water mass [kg]'\n";
+    gp << "set xlabel 'time[h]'\n";
     gp << "set ylabel ''\n";
   }
   else if (plt == "rwp")
   {
     gp << "set title 'rain water path [g / m^2]'\n";
-    gp << "set xlabel ''\n";
+    gp << "set xlabel 'time [h]'\n";
+    gp << "set ylabel ''\n";
+  }
+  else if (plt == "cwp")
+  {
+    gp << "set title 'cloud water path [g / m^2]'\n";
+    gp << "set xlabel 'time [h]'\n";
     gp << "set ylabel ''\n";
   }
   else if (plt == "cloud_base_dycoms")
diff --git a/drawbicyc/include/plot_fields.hpp b/drawbicyc/include/plot_fields.hpp
index 81f66c6..e53738d 100644
--- a/drawbicyc/include/plot_fields.hpp
+++ b/drawbicyc/include/plot_fields.hpp
@@ -19,6 +19,9 @@ void plot_fields(Plotter_t plotter, Plots plots, std::string type)
 
   auto& n = plotter.map;
 
+  // read in density
+  typename Plotter_t::arr_t rhod(plotter.h5load(plotter.file + "/const.h5", "G"));
+  
   blitz::firstIndex i;
   blitz::secondIndex j;
   blitz::Range all = blitz::Range::all();
@@ -71,6 +74,63 @@ void plot_fields(Plotter_t plotter, Plots plots, std::string type)
         }
         catch(...){}
       }
+      else if (plt == "mass_conc_cloud_droplets")
+      {
+	try{
+	// mass concentration of cloud droplets qc
+	auto tmp = plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"]) * 4. / 3 * 3.14 * 1e3 / rhod;
+	std::string title = " mass concentration of cloud droplets q_c";
+        gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+        plotter.plot(gp, tmp);
+        }
+        catch(...){}
+      }
+      else if (plt == "mass_conc_rain_drops")
+      {
+        try{
+        // mass concentration of rain droplets qr
+        auto tmp = plotter.h5load_timestep("rain_rw_mom3", at * n["outfreq"]) * 4. / 3 * 3.14 * 1e3 / rhod;
+        std::string title = " mass concentration of rain drops q_r";
+        gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+        plotter.plot(gp, tmp);
+        }
+        catch(...){}
+      }
+      else if (plt == "num_conc_droplets")
+      {
+	try{
+        // cloud particle concentration
+        auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]) * 1e-6;
+        std::string title ="activated droplets spec. conc. [mg^{-1}]";
+        gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+	plotter.plot(gp, tmp);
+        }
+        catch(...){}
+      }
+      else if (plt == "num_conc_rain_drops")
+      {
+        try{
+	// rain particle concentration
+	auto tmp = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]) * 1e-6;
+        std::string title = "rain (r > 25um) drops spec. conc. [mg^{-1}]";
+        gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+	gp << "set logscale cb\n";
+        plotter.plot(gp, tmp);
+        gp << "unset logscale cb\n";
+        }
+        catch(...){}
+      }	
+      else if (plt == "LWC")
+      {
+	try{
+	// liquid water content
+	auto tmp = (plotter.h5load_rc_timestep(at * n["outfreq"]) + plotter.h5load_rr_timestep(at * n["outfreq"])) * 1e3;//cloud + rain
+        std::string title = "liquid water content [g/kg]";
+        gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+        plotter.plot(gp, tmp);
+        }
+        catch(...){}
+      }
       else if (plt == "r_dry")
       {
         try{
@@ -113,6 +173,65 @@ void plot_fields(Plotter_t plotter, Plots plots, std::string type)
         }
         catch(...){}
       }
+      else if (plt == "effective_radius")
+      {
+        try{
+  	// effective radius
+  	typename Plotter_t::arr_t tmp2(plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]));
+        typename Plotter_t::arr_t tmp3(plotter.h5load_timestep("actrw_rw_mom3", at * n["outfreq"]));
+        typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+        auto tmp4 = iscloudy_rc_rico(snap);
+
+        auto tmp5 = where(tmp2 > 0, tmp3 / tmp2, 0);
+        auto ratio = tmp5 * tmp4*1e6;
+
+        std::string title = "act. droplet effective radius [μm]"; 
+        gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+        plotter.plot(gp, ratio);
+        }
+        catch(...){}
+      }
+      else if (plt == "ratio_mean_volume_r_to_eff_r_cubed")
+      {
+	try{
+	//Ratio of mean radius cubed k
+	typename Plotter_t::arr_t tmp(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]));
+        typename Plotter_t::arr_t tmp2(plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]));
+        typename Plotter_t::arr_t tmp3(plotter.h5load_timestep("actrw_rw_mom3", at * n["outfreq"]));
+        typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+        auto tmp4 = iscloudy_rc_rico(snap);
+
+        auto tmp5 = where(tmp > 0, tmp2 * tmp2 * tmp2 / tmp, 0);
+        auto tmp6 = where(tmp3 > 0, tmp5 / tmp3 / tmp3, 0);
+        auto ratio = tmp6 * tmp4;
+
+        std::string title = "Ratio of mean radius cubed k";
+        gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+        plotter.plot(gp, ratio);
+        }
+        catch(...){}
+      }
+      else if (plt == "std_dev_droplet_size_dist")
+      {
+	try{
+	//Stadndard deviation of the size distribution
+	typename Plotter_t::arr_t tmp1(plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]));
+        typename Plotter_t::arr_t tmp2(plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]));
+        typename Plotter_t::arr_t tmp0(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]));
+        typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+        auto tmp3 = iscloudy_rc_rico(snap);
+
+        auto tmp = where(tmp0 > 0, tmp2 / tmp0 - tmp1 / tmp0 * tmp1 / tmp0, 0.);
+        auto tmp4 = where(tmp < 0 , 0, tmp);
+        auto tmp5 = sqrt(tmp4);
+        auto tmp6 = tmp5 * tmp3 * 1e6;
+
+        std::string title = "Standard deviation of the size distribution [um]";
+        gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n";
+        plotter.plot(gp, tmp6);
+        }
+        catch(...){}
+      }	
 /*
       else if (plt == "na")
       {
diff --git a/drawbicyc/include/plot_prof.hpp b/drawbicyc/include/plot_prof.hpp
index 2e7d086..389eabd 100644
--- a/drawbicyc/include/plot_prof.hpp
+++ b/drawbicyc/include/plot_prof.hpp
@@ -42,7 +42,8 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool
   std::cerr << int(n["dt"] * n["outfreq"]+0.5) << std::endl;
   int first_timestep =  vm["prof_start"].as<int>() / int(n["dt"] * n["outfreq"]+0.5);
   int last_timestep =  vm["prof_end"].as<int>() / int(n["dt"] * n["outfreq"]+0.5);
-
+  //int start_time =  vm["prof_start"].as<int>(); 
+  //int end_time =  vm["prof_end"].as<int>(); 
   // some ugly constants
   const double p_1000 = 100000.;
   const double L = 2.5e6;
@@ -68,30 +69,113 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool
     blitz::secondIndex j;
     typename Plotter_t::arr_t res(rhod.shape());
     typename Plotter_t::arr_t res_tmp(rhod.shape());
+    //typename Plotter_t::arr_t res_tmp1(rhod.shape());
     typename Plotter_t::arr_t res_tmp2(rhod.shape());
+    //typename Plotter_t::arr_t res_tmp3(rhod.shape());
+    //typename Plotter_t::arr_t res_tmp4(rhod.shape());
+    //typename Plotter_t::arr_t res_sum(rhod.shape());
+    //typename Plotter_t::arr_t res_mean(rhod.shape());
+    //typename Plotter_t::arr_t rc_mask(rhod.shape());
     blitz::Array<float, 1> res_prof_sum(n["z"]);  // profile interpolate to the uniform grid summed over timesteps
     blitz::Array<float, 1> res_prof(n["z"]);      // profile interpolate to the uniform grid
     blitz::Array<float, 1> res_prof_hlpr(n["z"]); // actual profile
     blitz::Array<float, 1> prof_tmp(n["z"]);
+    //blitz::Array<float, 1> res_prof_num(n["z"]);
     blitz::Array<int, 1>   occur_no(n["z"]);      // number of occurances - for unusual profiles like base_prflux_vs_clhght
     blitz::Range all = blitz::Range::all();
 
     res_prof_sum = 0;
     occur_no = 0;
-
+    
+    if (plt == "coal_tele")
+      {
+	//coal tele mass flux in W/m^2
+	  //res = plotter.h5load_timestep("coal_tele_mass_flux",first_timestep *  n["outfreq"]) * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * first_timestep * n["outfreq"] * n["dt"]);
+	  res = plotter.h5load_timestep("coal_tele_mass_flux",last_timestep *  n["outfreq"]) * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * last_timestep * n["outfreq"] * n["dt"]);
+	  res -=  plotter.h5load_timestep("coal_tele_mass_flux",first_timestep *  n["outfreq"]) * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * first_timestep * n["outfreq"] * n["dt"]);
+	  // - plotter.h5load_timestep("coal_tele_mass_flux",first_timestep *  n["outfreq"]) * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * first_timestep * n["outfreq"] * n["dt"]);
+	  res_prof_hlpr = plotter.horizontal_mean(res);
+	  //auto tmp2 = plotter.h5load_timestep("coal_tele_mass_flux", last_timestep * n["outfreq"]);
+          //typename Plotter_t::arr_t snap(tmp);
+          //typename Plotter_t::arr_t snap2(tmp2);
+	  //res_prof_sum = (res_prof * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * last_timestep * n["outfreq"] * n["dt"])) - (res_prof_hlpr * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * first_timestep * n["outfreq"] * n["dt"]));  
+	  //res = tmp;
+      }
     for (int at = first_timestep; at <= last_timestep; ++at) // TODO: mark what time does it actually mean!
     {
       std::cout << at * n["outfreq"] << std::endl;
       res = 0;
 
+//      if (plt == "rliq")
+//      {
+	// liquid water content
+ //       res += plotter.h5load_ra_timestep(at * n["outfreq"]) * 1e3; // aerosol
+//        res += plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud
+//       res += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // rain
+//        res_prof_hlpr = plotter.horizontal_mean(res); // average in x
+//      }
       if (plt == "rliq")
       {
 	// liquid water content
- //       res += plotter.h5load_ra_timestep(at * n["outfreq"]) * 1e3; // aerosol
         res += plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud
         res += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // rain
-        res_prof_hlpr = plotter.horizontal_mean(res); // average in x
+        typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+        res_tmp = iscloudy_rc_rico(snap);
+        res *= res_tmp;
+        prof_tmp = plotter.horizontal_sum(res_tmp);
+        res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / prof_tmp, 0);
+      }
+      if (plt == "actrw_rw_cl_conc")
+      {
+        // concentration of activated droplets (r > rc) in cloudy cells
+        {
+          typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+          res_tmp = iscloudy_rc_rico(snap);
+        }
+        {
+          auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]);
+          typename Plotter_t::arr_t snap(tmp);
+          res_tmp2 = snap;
+          res_tmp2 *= rhod / 1e6; // per cm^3
+        }
+        // cloudy only
+        prof_tmp = plotter.horizontal_sum(res_tmp); 
+        res_tmp2 *= res_tmp;
+        res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp2) / prof_tmp, 0);
       }
+/*      if (plt == "cloud_water_std")
+      {
+        // cloud_water_std
+        //auto res = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud
+        //typename Plotter_t::arr_t snap(res);
+	res = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud
+	res_prof_sum = plotter.horizontal_sum(res);
+	res *=  iscloudy_rc_rico(res);
+	res_tmp = plotter.horizontal_sum(res_tmp) / res_prof_sum;
+	res_tmp2 = plotter.horizontal_sum((res - res_tmp) * (res - res_tmp)) / (res_prof_sum -1);
+	res_prof_hlpr = where(res >0, sqrt(res_tmp2), 0);	
+		
+	//rc_mask = iscloudy_rc_rico(res);
+        //res_tmp = res * rc_mask;
+        //res_prof_num = plotter.horizontal_sum(rc_mask);
+        //res_mean = plotter.horizontal_sum(res_tmp) / res_prof_num;
+        //res_sum = plotter.horizontal_sum((res_tmp - res_mean) * (res_tmp - res_mean)) / (res_prof_num -1);
+        //res_prof_hlpr = where(rc_mask >0 , sqrt(res_sum) , 0);
+      }
+*/
+      //if (plt == "rain_water_std")
+      //{
+        // rain_water_std
+      //  auto tmp = plotter.h5load_rr_timestep(at * n["outfreq"])
+
+      //  res = plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // cloud
+      //  rc_mask = iscloudy_rc_rico(res);
+      //  res_tmp = res * rc_mask;
+      //  res_prof_num = plotter.horizontal_sum(rc_mask);
+      //  res_mean = plotter.horizontal_sum(res_tmp) / res_prof_num;
+      //  res_sum = plotter.horizontal_sum((res_prof - res_mean) * (res_prof - res_mean)) / (res_prof_num -1);
+      //  res_prof_hlpr = where(rc_mask >0 , sqrt(res_sum) , 0);
+      //}
       if (plt == "gccn_conc")
       {
         res = plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"]) * rhod / 1e6;
@@ -420,7 +504,29 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool
         {
           auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]);
           typename Plotter_t::arr_t snap(tmp);
-          res_tmp = where(res_tmp > 0 , res_tmp / snap, res_tmp);
+          res_tmp = where(snap > 0 , res_tmp / snap, 0);
+        }
+        {
+          typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+          res_tmp2 = iscloudy_rc_rico(snap);
+          res_tmp *= res_tmp2;
+        }
+        // mean only over downdraught cells
+        prof_tmp = plotter.horizontal_sum(res_tmp2); // number of downdraft cells on a given level
+        res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp) / prof_tmp, 0);
+      }
+      if (plt == "actrw_rw_cl_sigma")
+      {
+        // sigma(radius) of activated droplets (r > rc) in cloudy cells
+        {
+          auto tmp = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]) * 1e6;
+          typename Plotter_t::arr_t snap(tmp);
+          res_tmp = snap; 
+        }
+        {
+          auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]);
+          typename Plotter_t::arr_t snap(tmp);
+          res_tmp = where(snap > 0 , res_tmp / snap, 0);
         }
         {
           typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
@@ -431,6 +537,63 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool
         prof_tmp = plotter.horizontal_sum(res_tmp2); // number of downdraft cells on a given level
         res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp) / prof_tmp, 0);
       }
+      //if (plt == "ratio_mean_volume_r_to_eff_r_cubed")
+      //{
+        // ratio mean (r volume / r effective) ^3 
+	//{
+	//   typename Plotter_t::arr_t snap(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]));
+        //    res_tmp = snap;
+        //  }
+        //  {
+        //    auto tmp2 = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]);
+        //    typename Plotter_t::arr_t snap(tmp2);
+        //    res_tmp2 = snap;
+        //   }
+        //   {
+        //    auto tmp3 = plotter.h5load_timestep("actrw_rw_mom3", at * n["outfreq"]);
+        //   typename Plotter_t::arr_t snap(tmp3);
+        //    res_tmp3 = snap;
+        //   }
+        //   {
+        //    typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+        //    res_tmp4 = iscloudy_rc_rico(snap);
+        //   }
+
+        //  res_tmp = where(res_tmp > 0, res_tmp2 * res_tmp2 * res_tmp2 / res_tmp, 0. );
+        //   res_tmp = where(res_tmp3 > 0, res_tmp / res_tmp3 / res_tmp3, 0.);
+        //   res_tmp *= res_tmp4;	 
+        //   prof_tmp = plotter.horizontal_sum(res_tmp4); // number of downdraft cells on a given level
+        //   res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp) / prof_tmp, 0);
+     // }
+     /* if (plt == "cloud_std_dev")
+      {
+          {
+            auto tmp1 = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]) *1e6;
+            typename Plotter_t::arr_t snap(tmp1);
+            res_tmp1 = snap;
+          }
+          {
+            auto tmp2 = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]) * 1e12;
+            typename Plotter_t::arr_t snap(tmp2);
+            res_tmp2 = snap;
+          }
+          {
+            auto tmp0 = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]);
+            typename Plotter_t::arr_t snap(tmp0);
+            res_tmp = snap;
+          }
+          {
+            typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+            res_tmp3 = iscloudy_rc_rico(snap);
+          }
+            res_tmp = where(res_tmp > 0, res_tmp2 / res_tmp - res_tmp1 / res_tmp * res_tmp1 / res_tmp, 0.);
+            res_tmp = where(res_tmp < 0 , 0, res_tmp);
+            res_tmp = sqrt(res_tmp);
+            res_tmp *= res_tmp3;
+            prof_tmp = plotter.horizontal_sum(res_tmp3); // number of downdraft cells on a given level
+            res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp) / prof_tmp, 0);
+      }
+	*/
       if (plt == "nc_up")
       {
         // updraft only
@@ -776,9 +939,14 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool
       {
 	// precipitation flux(doesnt include vertical volicty w!)
         { 
-          res = plotter.h5load_prflux_timestep(at * n["outfreq"]);
+          //res = plotter.dv;
+          //std::cerr << plotter.h5load_prflux_timestep(at * n["outfreq"]);
+	  //res = (plotter.dv).shape();
+          res = plotter.h5load_prflux_timestep(at * n["outfreq"]); // plotter.dv;
           res_prof_hlpr = plotter.horizontal_mean(res); // average in x
-        }
+	  //res_prof_hlpr = plotter.horizontal_sum(res)/122.; // average in x
+          
+	}
 	// add vertical velocity to precipitation flux (3rd mom of cloud drops * w)
 /*
         { 
@@ -929,11 +1097,95 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool
         res = where(res > 0 , res / res_tmp, res);
         res_prof_hlpr = plotter.horizontal_mean(res); // average in x
       }
+      if (plt == "actrw_mom1")
+      {
+        // mean radius of actrw droplets
+	res = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]) * 1e6;
+        res_tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]);
+        res = where(res > 0 , res / res_tmp, res);
+        res_prof_hlpr = plotter.horizontal_mean(res); // average in x
+      }
+      else if (plt == "actrw_mom2")
+      {
+        // variance of actrw droplet size distribution
+        res = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]) * 1e6 * 1e3; // convert to (um)^2 * 1e3
+	res_tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]);
+        res = where(res > 0 , res / res_tmp, res);
+        res_prof_hlpr = plotter.horizontal_mean(res); // average in x
+      }
+      else if (plt == "disp_r")
+      {
+        // relative dispersion (std dev / mean) of droplet radius distribution averaged over cells away from walls
+        //typename Plotter_t::arr_t m0(plotter.nowall(typename Plotter_t::arr_t(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])), distance_from_walls));
+        typename Plotter_t::arr_t m0(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]));
+        typename Plotter_t::arr_t m1(plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]));
+        typename Plotter_t::arr_t m2(plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]));
 
+        // calculate stddev of radius, store in m2
+        m2 = where(m0 > 0,
+            m2 / m0 - m1 / m0 * m1 / m0 , 0.);
 
+        // might be slightly negative due to numerical errors
+        m2 = where(m2 < 0, 0, m2);
+        m2 = sqrt(m2); // sqrt(variance)
 
-// ======================================================================================
+        //calculate mean radius, store in m1
+        m1 = where(m0 > 0,
+            m1 / m0, 0.);
+        
+        res = where(m1 > 0 , m2 / m1 , 0);
+
+        typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+        res_tmp = iscloudy_rc_rico(snap);
+        res *= res_tmp;
+        prof_tmp = plotter.horizontal_sum(res_tmp);
+
+        res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / prof_tmp, 0);
+     }
+     else if (plt == "sigma_r")
+     {
+        // std dev of droplet radius distribution averaged over cells 
+        //typename Plotter_t::arr_t m0(plotter.nowall(typename Plotter_t::arr_t(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])), distance_from_walls));
+        typename Plotter_t::arr_t m0(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]));
+        typename Plotter_t::arr_t m1(plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]));
+        typename Plotter_t::arr_t m2(plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]));
+
+        // calculate stddev of radius, store in m2
+        m2 = where(m0 > 0,
+            m2 / m0 - m1 / m0 * m1 / m0 , 0.);
 
+        // might be slightly negative due to numerical errors
+        m2 = where(m2 < 0, 0, m2);
+        res = sqrt(m2); // sqrt(variance)
+        
+        typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+        res_tmp = iscloudy_rc_rico(snap);
+        res *= res_tmp;
+        prof_tmp = plotter.horizontal_sum(res_tmp);
+
+        res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / prof_tmp, 0);
+     }
+     else if (plt == "mean_r")
+     {
+        // relative dispersion (std dev / mean) of droplet radius distribution averaged over cells away from walls
+        //typename Plotter_t::arr_t m0(plotter.nowall(typename Plotter_t::arr_t(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])), distance_from_walls));
+        typename Plotter_t::arr_t m0(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]));
+        typename Plotter_t::arr_t m1(plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]));
+
+        //calculate mean radius, store in m1
+        res = where(m0 > 0,
+            m1 / m0, 0.);
+        
+        typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+        res_tmp = iscloudy_rc_rico(snap);
+        res *= res_tmp;
+        prof_tmp = plotter.horizontal_sum(res_tmp);
+
+        res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / prof_tmp, 0);
+     }
+
+// ======================================================================================
+	
       if(normalize)
       {
         // interpolate profile to uniform vertical grid (height normalized by inversion height)
@@ -975,7 +1227,6 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool
       oprof_file << res_pos;
       res_pos_out_done = true;
     }
-
     if (plt != "base_prflux_vs_clhght")
       res_prof_sum /= last_timestep - first_timestep + 1;
     else
@@ -987,7 +1238,6 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool
     // do the plotting
     gp << "plot '-' with line\n";
     gp.send1d(boost::make_tuple(res_prof_sum, res_pos));
-
     if (plt == "base_prflux_vs_clhght")
     {
       oprof_file << plt << " number of occurances" << endl;
diff --git a/drawbicyc/include/plot_series.hpp b/drawbicyc/include/plot_series.hpp
index d2ac874..b1f536e 100644
--- a/drawbicyc/include/plot_series.hpp
+++ b/drawbicyc/include/plot_series.hpp
@@ -28,7 +28,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
   typename Plotter_t::arr_t rtot(rhod.shape());
 
   typename Plotter_t::arr_t res_tmp(rhod.shape());
-
+  
   // for calculating running averages of u and w, needed in TKE calc in Pi chamber LES
   std::vector<typename Plotter_t::arr_t> prev_u_vec, prev_w_vec;
   // container for the running sum
@@ -176,6 +176,27 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
       }
+      // cloud top height
+      else if (plt =="cloud_top_height")
+      {
+        try
+        {
+	  // cloud fraction (cloudy if ql > 1e-5))
+	  auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]);
+          typename Plotter_t::arr_t snap(tmp);
+          snap = iscloudy_rc_rico(snap);
+          plotter.k_i = blitz::last((snap == 1), plotter.LastIndex);
+          auto cloudy_column = plotter.k_i.copy();
+          cloudy_column = blitz::sum(snap, plotter.LastIndex);
+          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) = blitz::max(plotter.k_i)*n["dz"];
+         else
+           res_series[plt](at) = 0;
+        }
+        catch(...) {if(at==first_timestep) data_found[plt]=0;}
+      }
       // average sd_conc
       else if (plt == "sd_conc")
       {
@@ -512,6 +533,17 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
       }
+      else if (plt == "total_droplets_number")
+      {
+        // Total number of Cloud and Rain Droplets
+        try
+        {
+          auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at *n["outfreq"]) * rhod;
+	  typename Plotter_t::arr_t snap(tmp);
+	  res_series[plt](at) = blitz::sum(snap)*plotter.CellVol;
+	}
+      catch(...) {if(at==first_timestep) data_found[plt]=0;} 
+      }
       else if (plt == "cloud_base_dycoms")
       {
         // average cloud base in the domain
@@ -667,6 +699,15 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
       }
+      else if (plt == "acc_vol_precip")
+      {
+        // accumulated surface precipitation [m^3]
+        try
+        {
+          res_series[plt](at) = plotter.calc_acc_surf_precip_volume(prec_vol);
+        }
+        catch(...) {if(at==first_timestep) data_found[plt]=0;}
+      }
       else if (plt == "cl_acnv25_dycoms")
       {
         // autconversion rate with rain threshold r=25um (cloudy cells as in Dycoms) [g/(m3*s)]
@@ -794,7 +835,69 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
           }
         }
         catch(...) {if(at==first_timestep) data_found[plt]=0;}
-      }   
+      }       
+      else if (plt == "cwp")
+      {   
+        // cloud water path
+        try
+        {
+          {
+            typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"]));
+            snap *= rhod * 1e3; // water per cubic metre (should be wet density...) & g/kg
+            res_series[plt](at) = blitz::mean(snap); 
+          }
+        }
+        catch(...) {if(at==first_timestep) data_found[plt]=0;}
+      }       
+      else if (plt == "lwm")
+      {
+        //liquid water mass
+        try
+        {
+          { 
+            auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * rhod;
+            typename Plotter_t::arr_t snap(tmp);
+            snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * rhod;
+            snap *= plotter.CellVol;
+            snap(plotter.hrzntl_slice(0)) = snap(plotter.hrzntl_slice(0))/2;
+            snap(plotter.hrzntl_slice(-1)) = snap(plotter.hrzntl_slice(-1))/2;
+            res_series[plt](at) = blitz::sum(snap);
+          }
+        }
+        catch(...) {if(at==first_timestep) data_found[plt]=0;}
+      }
+      else if (plt == "cwm")
+      {
+	//cloud water mass
+	try
+	{
+	  { 
+            auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]);
+            typename Plotter_t::arr_t snap(tmp);
+            snap *= rhod * plotter.CellVol;
+            snap(plotter.hrzntl_slice(0)) = snap(plotter.hrzntl_slice(0))/2;
+            snap(plotter.hrzntl_slice(-1)) = snap(plotter.hrzntl_slice(-1))/2;
+            res_series[plt](at) = blitz::sum(snap);
+          }
+	}
+	catch(...) {if(at==first_timestep) data_found[plt]=0;}
+      }
+      else if (plt == "rwm")
+      {
+	//liquid water mass
+	try
+	{
+	  {
+	    auto tmp = plotter.h5load_rr_timestep(at * n["outfreq"]);
+            typename Plotter_t::arr_t snap(tmp);
+            snap *= rhod * plotter.CellVol;
+            snap(plotter.hrzntl_slice(0)) = snap(plotter.hrzntl_slice(0))/2;
+            snap(plotter.hrzntl_slice(-1)) = snap(plotter.hrzntl_slice(-1))/2;
+            res_series[plt](at) = blitz::sum(snap);
+	  }
+	}
+	catch(...) {if(at==first_timestep) data_found[plt]=0;}
+      }
       else if (plt == "surf_flux_latent")
       {   
         try
@@ -1197,15 +1300,16 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
       else if (plt == "cl_meanr")
       {
         // cloud droplets mean radius in cloudy grid cells
+        // modified for actrw_mom1. cloud -> actrw
         try
         {
           // cloud fraction (cloudy if N_c > 20/cm^3)
-          typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]));
+          typename Plotter_t::arr_t snap(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]));
           snap *= rhod; // b4 it was specific moment
           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"]));
-          typename Plotter_t::arr_t snap_m1(plotter.h5load_timestep("cloud_rw_mom1", at * n["outfreq"])*1e6); // in microns
+          typename Plotter_t::arr_t snap_m0(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]));
+          typename Plotter_t::arr_t snap_m1(plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"])*1e6); // in microns
           snap_m0 *= snap;
           snap_m1 *= snap;
           auto tot_m0 = blitz::sum(snap_m0);
@@ -1773,6 +1877,10 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
     {
       res_pos *= 60.;
     }
+    else if (plt == "cloud_top_height")
+    {
+      res_pos *= 60.;
+    }
     else if (plt == "tot_tke" || plt == "tot_tke" || plt == "sgs_tke" || plt == "uw_resolved_tke")
     {
       res_pos *= 60.;
@@ -1785,6 +1893,10 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type)
     {
       res_series[plt] *= (n["z"] - 1) * n["dz"]; // top and bottom cells are smaller
     }
+    else if (plt == "cwp")
+    {
+      res_series[plt] *= (n["z"] - 1) * n["dz"]; // top and bottom cells are smaller
+    }
     else if (plt == "er")
     {
       // central difference, in cm
diff --git a/drawbicyc/include/plots.hpp b/drawbicyc/include/plots.hpp
index 619b592..db99a69 100644
--- a/drawbicyc/include/plots.hpp
+++ b/drawbicyc/include/plots.hpp
@@ -5,9 +5,10 @@
 #include "cases/moist_thermal/plots.hpp"
 #include "cases/PiChamber/plots.hpp"
 #include "cases/PiChamber/plots_icmw.hpp"
-#include "cases/Lasher_Trapp/plots.hpp"
+#include "cases/Cumulus_Congestus/plots.hpp"
 #include "cases/common_plots/GCCN_CCN_conc/plots.hpp"
 
+
 const std::vector<std::string> series_sgs({
 // "tot_tke"
 });
@@ -53,10 +54,10 @@ class Plots
       series.insert(series.end(), series_PiChamber.begin(), series_PiChamber.end());
       fields.insert(fields.end(), fields_PiChamber.begin(), fields_PiChamber.end());
     }
-    else if(type == "Lasher_Trapp") {
-     // profs.insert(profs.end(), profs_Lasher_Trapp.begin(), profs_Lasher_Trapp.end());
-      series.insert(series.end(), series_Lasher_Trapp.begin(), series_Lasher_Trapp.end());
-      fields.insert(fields.end(), fields_Lasher_Trapp.begin(), fields_Lasher_Trapp.end());
+    else if(type == "cumulus_congestus") {
+      profs.insert(profs.end(), profs_Cumulus_Congestus.begin(), profs_Cumulus_Congestus.end());
+      series.insert(series.end(), series_Cumulus_Congestus.begin(), series_Cumulus_Congestus.end());
+      fields.insert(fields.end(), fields_Cumulus_Congestus.begin(), fields_Cumulus_Congestus.end());
     }
     else if(type == "pi_chamber_icmw") { 
       profs.insert(profs.end(), profs_PiChamberICMW.begin(), profs_PiChamberICMW.end());