diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp
index 949f94ab..abb41acf 100644
--- a/bindings/python/lib.cpp
+++ b/bindings/python/lib.cpp
@@ -131,10 +131,13 @@ BOOST_PYTHON_MODULE(libcloudphxx)
     bp::def("exner", &common::exner<real_t>);
     bp::def("p_v", &common::p_v<real_t>);
     bp::def("p_vs", &common::p_vs<real_t>);
+    bp::def("p_vsi", &common::p_vsi<real_t>);
     bp::def("r_vs", &common::r_vs<real_t>);
     bp::def("r_vsi", &common::r_vsi<real_t>);
     bp::def("p_vs_tet", &common::p_vs_tet<real_t>);
     bp::def("l_v", &common::l_v<real_t>);
+    bp::def("l_s", &common::l_s<real_t>);
+    bp::def("l_f", &common::l_f<real_t>);
     bp::def("T", &common::T<real_t>);
     bp::def("p", &common::p<real_t>);
     bp::def("visc", &common::visc<real_t>);
@@ -282,6 +285,7 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("turb_cond", &lgr::opts_t<real_t>::turb_cond)
       .def_readwrite("turb_coal", &lgr::opts_t<real_t>::turb_coal)
       .def_readwrite("ice_nucl", &lgr::opts_t<real_t>::ice_nucl)
+      .def_readwrite("depo", &lgr::opts_t<real_t>::depo)
       .def_readwrite("dt", &lgr::opts_t<real_t>::dt)
       .add_property("src_dry_distros", &lgrngn::get_sdd<real_t>, &lgrngn::set_sdd<real_t>)
       .add_property("src_dry_sizes", &lgrngn::get_ds<real_t>, &lgrngn::set_sds<real_t>)
@@ -427,7 +431,12 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def("diag_water_cons",    &lgr::particles_proto_t<real_t>::diag_water_cons)
       .def("diag_ice_a_mom",    &lgr::particles_proto_t<real_t>::diag_ice_a_mom)
       .def("diag_ice_c_mom",    &lgr::particles_proto_t<real_t>::diag_ice_c_mom)
+      .def("diag_ice_a_rng",    &lgr::particles_proto_t<real_t>::diag_ice_a_rng)
+      .def("diag_ice_c_rng",    &lgr::particles_proto_t<real_t>::diag_ice_c_rng)
+      .def("diag_ice_a_rng_cons",    &lgr::particles_proto_t<real_t>::diag_ice_a_rng_cons)
+      .def("diag_ice_c_rng_cons",    &lgr::particles_proto_t<real_t>::diag_ice_c_rng_cons)
       .def("diag_ice_mix_ratio",    &lgr::particles_proto_t<real_t>::diag_ice_mix_ratio)
+      .def("diag_sstp_cond_mom",    &lgr::particles_proto_t<real_t>::diag_sstp_cond_mom)
       .def("outbuf",       &lgrngn::outbuf<real_t>)
       .def("get_attr",    &lgr::particles_proto_t<real_t>::get_attr)
     ;
diff --git a/include/libcloudph++/lgrngn/particles.hpp b/include/libcloudph++/lgrngn/particles.hpp
index c0546e54..6157e067 100644
--- a/include/libcloudph++/lgrngn/particles.hpp
+++ b/include/libcloudph++/lgrngn/particles.hpp
@@ -118,6 +118,7 @@ namespace libcloudphxx
       virtual void diag_vp_mom(const int&)                                      { assert(false); }
       virtual void diag_wp_mom(const int&)                                      { assert(false); }
       virtual void diag_incloud_time_mom(const int&)                            { assert(false); } // requires opts_init.diag_incloud_time==true
+      virtual void diag_sstp_cond_mom(const int &k)                             { assert(false); }
       virtual void diag_max_rw()                                                { assert(false); }
       virtual void diag_vel_div()                                               { assert(false); }
       virtual std::map<libcloudphxx::common::output_t, real_t> diag_puddle()    { assert(false); return std::map<libcloudphxx::common::output_t, real_t>(); }
@@ -212,6 +213,7 @@ namespace libcloudphxx
       void diag_vp_mom(const int&);
       void diag_wp_mom(const int&);
       void diag_incloud_time_mom(const int &k);
+      void diag_sstp_cond_mom(const int &k);
       void diag_wet_mass_dens(const real_t&, const real_t&);
 
       void diag_chem(const enum common::chem::chem_species_t&);
@@ -321,6 +323,7 @@ namespace libcloudphxx
       void diag_vp_mom(const int&);
       void diag_wp_mom(const int&);
       void diag_incloud_time_mom(const int&);
+      void diag_sstp_cond_mom(const int &k);
       void diag_wet_mass_dens(const real_t&, const real_t&);
       std::vector<real_t> get_attr(const std::string &);
       real_t *outbuf();
diff --git a/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp b/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp
index 9dbf2401..097ddb50 100644
--- a/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp
+++ b/src/impl/condensation_deposition/perparticle/release_arrays_for_perparticle_sstp.ipp
@@ -30,7 +30,7 @@ namespace libcloudphxx
 
       if(opts_init.adaptive_sstp_cond)
       {
-        perparticle_sstp_cond_gp.reset();
+        // perparticle_sstp_cond_gp.reset();
       }
     }
   };
diff --git a/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp b/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp
index aeff9f63..b811501a 100644
--- a/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp
+++ b/src/impl/diagnose_SD_attributes/particles_impl_moms.ipp
@@ -396,5 +396,57 @@ namespace libcloudphxx
     {
       moms_calc(vec_bgn, n_part, power, specific);
     }
+
+    // like moms_calc, but not accounting for multiplicity (n) and not dividing by volume (or air mass)
+    template <typename real_t, backend_t device>
+    template <typename it_t> // iterator type
+    void particles_t<real_t, device>::impl::SD_moms_calc( 
+      const it_t &vec_bgn,
+      const thrust_size_t npart,
+      const real_t power
+    )
+    {
+      thrust::pair<
+        thrust_device::vector<thrust_size_t>::iterator,
+        typename thrust_device::vector<real_t>::iterator
+      > it_pair = thrust::reduce_by_key(
+        // input - keys
+        sorted_ijk.begin(), sorted_ijk.begin()+npart,  
+        // input - values
+        thrust::make_transform_iterator(
+          thrust::make_zip_iterator(thrust::make_tuple(
+            thrust::make_constant_iterator<real_t>(1),
+            thrust::make_permutation_iterator(vec_bgn, sorted_id.begin())
+          )),
+          detail::moment_counter<real_t>(power)
+        ),
+        // output - keys
+        count_ijk.begin(),
+        // output - values
+        count_mom.begin()
+      );  
+
+      count_n = it_pair.first - count_ijk.begin();
+#if !defined(NDEBUG)
+      {
+        int nan_count = thrust::transform_reduce(count_mom.begin(), count_mom.begin() + count_n, isnaninf(), 0, thrust::plus<bool>());
+        if(nan_count>0)
+        {
+          std::cout << nan_count << " nan/inf numbers detected in count_mom after reduce_by_key " << std::endl;
+        }
+      }
+#endif
+      assert(count_n <= n_cell);
+    }
+
+    template <typename real_t, backend_t device>
+    template <typename it_t> // iterator type
+    void particles_t<real_t, device>::impl::SD_moms_calc(
+      const it_t &vec_bgn,
+      const real_t power
+    )
+    {
+      moms_calc(vec_bgn, n_part, power);
+    }
   };  
 };
diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp
index f4c9e8f0..c254b802 100644
--- a/src/impl/particles_impl.ipp
+++ b/src/impl/particles_impl.ipp
@@ -516,7 +516,7 @@ namespace libcloudphxx
         if(allow_sstp_cond && opts_init.exact_sstp_cond && sstp_cond_exact_nomix_adaptive)
           tmp_drp_no = std::max(tmp_drp_no, 4); // why 5? not 4?
         if(allow_sstp_cond && opts_init.exact_sstp_cond && !sstp_cond_exact_nomix_adaptive)
-          tmp_drp_no = std::max(tmp_drp_no, 7); // why 8? not 7?
+          tmp_drp_no = std::max(tmp_drp_no, 10); // for some reason it fails for less than 10
         // if(allow_sstp_cond && opts_init.exact_sstp_cond && opts_init.const_p)
         //   tmp_drp_no = std::max(tmp_drp_no, 7);
         tmp_device_real_part.add_vectors(tmp_drp_no-1); // -1 because 1 is already created in the ctor
@@ -672,6 +672,17 @@ namespace libcloudphxx
         const real_t power,
         const bool specific = true
       );
+      template<typename it_t> // iterator type
+      void SD_moms_calc(
+        const it_t &vec_bgn,
+        const thrust_size_t npart,
+        const real_t power
+      );
+      template<typename it_t> // iterator type
+      void SD_moms_calc(
+        const it_t &vec_bgn,
+        const real_t power
+      );
 
       void mass_dens_estim(
         const typename thrust_device::vector<real_t>::iterator &vec_bgn,
diff --git a/src/particles_diag.ipp b/src/particles_diag.ipp
index 720a3935..1c7c1e36 100644
--- a/src/particles_diag.ipp
+++ b/src/particles_diag.ipp
@@ -654,5 +654,19 @@ namespace libcloudphxx
     {
       return pimpl->output_puddle;
     }
+
+    template <typename real_t, backend_t device>
+    void particles_t<real_t, device>::diag_sstp_cond_mom(const int &n)
+    {       
+      if(pimpl->opts_init.exact_sstp_cond && (pimpl->sstp_cond > 1 || pimpl->sstp_cond_act > 1) && pimpl->opts_init.adaptive_sstp_cond) 
+      {
+        if (pimpl->perparticle_sstp_cond_gp)
+        {
+          pimpl->SD_moms_calc(pimpl->perparticle_sstp_cond_gp->get().begin(), n);
+        }
+      }
+      else
+        assert(0 && "diag_sstp_cond_mom called, but adaptive substepping is off (opts_init.exact_sstp_cond && (opts_ini.sstp_cond > 1 || opts_ini.sstp_cond_act > 1) && opts_init.adaptive_sstp_cond) == False. Therefore number of substeps is defined by opts_init.sstp_cond.");
+    }   
   };
 };
diff --git a/src/particles_multi_gpu_diag.ipp b/src/particles_multi_gpu_diag.ipp
index 797b134f..2f12a33f 100644
--- a/src/particles_multi_gpu_diag.ipp
+++ b/src/particles_multi_gpu_diag.ipp
@@ -213,6 +213,12 @@ namespace libcloudphxx
       pimpl->mcuda_run(&particles_t<real_t, CUDA>::diag_incloud_time_mom, k);
     }   
 
+    template <typename real_t>
+    void particles_t<real_t, multi_CUDA>::diag_sstp_cond_mom(const int &k)
+    {
+      pimpl->mcuda_run(&particles_t<real_t, CUDA>::diag_sstp_cond_mom, k);
+    }
+    
     template <typename real_t>
     void particles_t<real_t, multi_CUDA>::diag_wet_mass_dens(const real_t &a, const real_t &b)
     {