diff --git a/bindings/python/blk_1m.hpp b/bindings/python/blk_1m.hpp
index 1d918900..db1c227f 100644
--- a/bindings/python/blk_1m.hpp
+++ b/bindings/python/blk_1m.hpp
@@ -121,48 +121,47 @@ namespace libcloudphxx
 			np2bz<arr_t>(rr)
 		  );
       }
-      
-      template <class arr_t>
-      void rhs_cellwise_ice(
-	const b1m::opts_t<typename arr_t::T_numtype> &opts,
-	bp_array &dot_rc,
-	bp_array &dot_rr,
-	bp_array &dot_rv,
-	bp_array &dot_ria,
-	const bp_array &rc,
-	const bp_array &rr,
-	const bp_array &rv,
-	const bp_array &ria,
-	const bp_array &theta,
-	const bp_array &p,
-	const bp_array &rhod,
-	const typename arr_t::T_numtype &dt
-      ) 
+
+    	template <class arr_t>
+		 void rhs_cellwise_nwtrph(
+	   const b1m::opts_t<typename arr_t::T_numtype> &opts,
+	   bp_array &dot_th,
+	   bp_array &dot_rv,
+	   bp_array &dot_rc,
+	   bp_array &dot_rr,
+	   const bp_array &rhod,
+	   const bp_array &p,
+	   const bp_array &th,
+	   const bp_array &rv,
+	   const bp_array &rc,
+	   const bp_array &rr,
+	   const typename arr_t::T_numtype &dt
+		 )
       {
-	arr_t
-	  np2bz_dot_rc(np2bz<arr_t>(dot_rc)), 
-	  np2bz_dot_rr(np2bz<arr_t>(dot_rr)),
-      np2bz_dot_rv(np2bz<arr_t>(dot_rv)),
-      np2bz_dot_ria(np2bz<arr_t>(dot_ria));
-	b1m::rhs_cellwise_ice(
-	  opts,
-	  np2bz_dot_rc,
-	  np2bz_dot_rr,
-	  np2bz_dot_rv,
-	  np2bz_dot_ria,
-	  np2bz<arr_t>(rc),
-	  np2bz<arr_t>(rr),
-	  np2bz<arr_t>(rv),
-	  np2bz<arr_t>(ria),
-	  np2bz<arr_t>(theta),
-	  np2bz<arr_t>(p),
-	  np2bz<arr_t>(rhod),
-	  dt
-	);
-      } 
+      	arr_t
+			np2bz_dot_rc(np2bz<arr_t>(dot_rc)),
+			np2bz_dot_rr(np2bz<arr_t>(dot_rr)),
+			np2bz_dot_rv(np2bz<arr_t>(dot_rv)),
+			np2bz_dot_th(np2bz<arr_t>(dot_th));
+      	b1m::rhs_cellwise_nwtrph(
+			opts,
+			np2bz_dot_th,
+			np2bz_dot_rv,
+			np2bz_dot_rc,
+			np2bz_dot_rr,
+			np2bz<arr_t>(rhod),
+			np2bz<arr_t>(p),
+			np2bz<arr_t>(th),
+			np2bz<arr_t>(rv),
+			np2bz<arr_t>(rc),
+			np2bz<arr_t>(rr),
+				dt
+		  );
+      }
+
 
       template <class arr_t>
-      void rhs_cellwise_nwtrph(
+      void rhs_cellwise_nwtrph_ice(
 	const b1m::opts_t<typename arr_t::T_numtype> &opts,
 	bp_array &dot_th,
 	bp_array &dot_rv,
@@ -185,7 +184,7 @@ namespace libcloudphxx
 	  np2bz_dot_rr(np2bz<arr_t>(dot_rr)),
 	  np2bz_dot_rv(np2bz<arr_t>(dot_rv)),
 	  np2bz_dot_th(np2bz<arr_t>(dot_th));
-	b1m::rhs_cellwise_nwtrph(
+	b1m::rhs_cellwise_nwtrph_ice(
 	  opts,
 	  np2bz_dot_th,
 	  np2bz_dot_rv,
diff --git a/bindings/python/lib.cpp b/bindings/python/lib.cpp
index d7adc5f4..dfe3b84c 100644
--- a/bindings/python/lib.cpp
+++ b/bindings/python/lib.cpp
@@ -155,7 +155,6 @@ BOOST_PYTHON_MODULE(libcloudphxx)
       .def_readwrite("conv", &b1m::opts_t<real_t>::conv)
       .def_readwrite("accr", &b1m::opts_t<real_t>::accr)
       .def_readwrite("sedi", &b1m::opts_t<real_t>::sedi)
-      .def_readwrite("ice", &b1m::opts_t<real_t>::ice)
       .def_readwrite("homA1", &b1m::opts_t<real_t>::homA1)
       .def_readwrite("homA2", &b1m::opts_t<real_t>::homA2)
       .def_readwrite("hetA", &b1m::opts_t<real_t>::hetA)
@@ -167,8 +166,8 @@ BOOST_PYTHON_MODULE(libcloudphxx)
     bp::def("adj_cellwise_constp", blk_1m::adj_cellwise_constp<arr_t>);
     bp::def("adj_cellwise_nwtrph", blk_1m::adj_cellwise_nwtrph<arr_t>);
     bp::def("rhs_cellwise", blk_1m::rhs_cellwise<arr_t>);
-    bp::def("rhs_cellwise_ice", blk_1m::rhs_cellwise_ice<arr_t>);
-    bp::def("rhs_cellwise_nwtrph", blk_1m::rhs_cellwise_nwtrph<arr_t>); 
+    bp::def("rhs_cellwise_nwtrph", blk_1m::rhs_cellwise_nwtrph<arr_t>);
+    bp::def("rhs_cellwise_nwtrph_ice", blk_1m::rhs_cellwise_nwtrph_ice<arr_t>);
     bp::def("rhs_columnwise", blk_1m::rhs_columnwise<arr_t>); // TODO: handle the returned flux
   }
 
diff --git a/include/libcloudph++/blk_1m/options.hpp b/include/libcloudph++/blk_1m/options.hpp
index 38526171..792b7e44 100644
--- a/include/libcloudph++/blk_1m/options.hpp
+++ b/include/libcloudph++/blk_1m/options.hpp
@@ -21,7 +21,6 @@ namespace libcloudphxx
         conv = true,    // autoconversion
         accr = true,    // accretion
         sedi = true,    // sedimentation
-        ice = true,    // enable ice processes
         homA1 = true,    // homogeneous nucleation 1 of ice A
         homA2 = true,    // homogeneous nucleation 2 of ice A
         hetA = true;    // heterogeneous nucleation of ice A
diff --git a/include/libcloudph++/blk_1m/rhs_cellwise.hpp b/include/libcloudph++/blk_1m/rhs_cellwise.hpp
index 1daf00e2..af5b2fa6 100644
--- a/include/libcloudph++/blk_1m/rhs_cellwise.hpp
+++ b/include/libcloudph++/blk_1m/rhs_cellwise.hpp
@@ -63,105 +63,72 @@ namespace libcloudphxx
       }
     }
 
-//<listing>
     template <typename real_t, class cont_t>
-    void rhs_cellwise_ice(
+    void rhs_cellwise_nwtrph(
       const opts_t<real_t> &opts,
+      cont_t &dot_th_cont,
+      cont_t &dot_rv_cont,
       cont_t &dot_rc_cont,
       cont_t &dot_rr_cont,
-      cont_t &dot_rv_cont,
-      cont_t &dot_ria_cont,
+      const cont_t &rhod_cont,
+      const cont_t &p_cont,
+      const cont_t &th_cont,
+      const cont_t &rv_cont,
       const cont_t &rc_cont,
       const cont_t &rr_cont,
-      const cont_t &rv_cont,
-      const cont_t &ria_cont,
-      const cont_t &th_cont,
-      const cont_t &p_cont,
-      const cont_t &rhod_cont,
       const real_t &dt
     )
-//</listing>
     {
-      for (auto tup : zip(dot_rc_cont, dot_rr_cont, dot_rv_cont, dot_ria_cont, rc_cont, rr_cont, rv_cont, ria_cont, th_cont, p_cont, rhod_cont))
+      rhs_cellwise<real_t, cont_t>(opts, dot_rc_cont, dot_rr_cont,  rc_cont, rr_cont);
+
+      // rain evaporation treated as a force in Newthon-Raphson saturation adjustment
+      for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rr_cont, rhod_cont, p_cont, th_cont, rv_cont, rr_cont))
       {
         using namespace common;
+
         real_t
-          rc_to_ria = 0,
-          rv_to_ria = 0,
-          &dot_rc = std::get<0>(tup),
-          &dot_rr = std::get<1>(tup),
-          &dot_rv = std::get<2>(tup),
-          &dot_ria = std::get<3>(tup);
-        const real_t
-          &rc     = std::get<4>(tup),
-          &rr     = std::get<5>(tup),
-          &rv     = std::get<6>(tup),
-          &ria     = std::get<7>(tup);
+          rr_to_rv = 0,
+          &dot_th = std::get<0>(tup),
+          &dot_rv = std::get<1>(tup),
+          &dot_rr = std::get<2>(tup);
 
-        const quantity<si::temperature, real_t>
-        th = std::get<8>(tup) * si::kelvins;
+        const quantity<si::mass_density, real_t>
+          rhod = std::get<3>(tup) * si::kilograms / si::cubic_metres;
 
         const quantity<si::pressure, real_t>
-        p   = std::get<9>(tup) * si::pascals;
-
-        const quantity<si::mass_density, real_t>
-        rhod = std::get<10>(tup) * si::kilograms / si::cubic_metres;
+          p   = std::get<4>(tup) * si::pascals;
 
-        quantity<si::temperature, real_t> T = th * theta_std::exner(p);
-        quantity<si::dimensionless, real_t> rvs = const_cp::r_vs(T, p);
-        quantity<si::dimensionless, real_t> rvsi = const_cp::r_vsi(T, p);
+        const quantity<si::temperature, real_t>
+          th = std::get<5>(tup) * si::kelvins;
 
-        // ice A heterogeneous nucleation
-        if (opts.hetA)
-        {
-          rc_to_ria += (
-            formulae::het_A_nucleation(
-                    ria * si::dimensionless(),
-                    rc * si::dimensionless(),
-                    th,
-                	p,
-                    rhod,
-                    dt * si::seconds
-                  ) * si::seconds // to make it dimensionless
-          );
-        }
+        const real_t
+          &rv = std::get<6>(tup),
+          &rr = std::get<7>(tup);
 
-        // ice A homogeneous nucleation 1
-        if (opts.homA1)
-        {
-         rv_to_ria += (
-           formulae::hom_A_nucleation_1(
-                   rv * si::dimensionless(),
-                   rvs * si::dimensionless(),
-                   rvsi * si::dimensionless(),
-                   th,
-                   p,
-                   dt * si::seconds
-                 ) * si::seconds // to make it dimensionless
-         );
-        }
+        quantity<si::temperature, real_t> T = th * theta_std::exner(p);
+        real_t r_vs = const_cp::r_vs(T, p);
 
-        // ice A homogeneous nucleation 2
-        if (opts.homA2)
-        {
-        rc_to_ria += (
-          formulae::hom_A_nucleation_2(
-                  rc * si::dimensionless(),
-                  th,
-                  p,
-                  dt * si::seconds
-                ) * si::seconds // to make it dimensionless
+        rr_to_rv += (
+          formulae::evaporation_rate(
+            rv * si::dimensionless(),
+            r_vs * si::dimensionless(),
+            rr * si::dimensionless(),
+            rhod,
+            p
+          ) * si::seconds * dt
         );
-        }
 
-	dot_rc -= rc_to_ria;
-    dot_rv -= rv_to_ria;
-    dot_ria += rc_to_ria + rv_to_ria;
+        // limiting
+        rr_to_rv = std::min(rr, rr_to_rv) / dt;
+
+        dot_rv += rr_to_rv;
+        dot_rr -= rr_to_rv;
+        dot_th -= const_cp::l_v(T) / (moist_air::c_pd<real_t>() * theta_std::exner(p)) * dot_rv / si::kelvins;
       }
     }
 
     template <typename real_t, class cont_t>
-    void rhs_cellwise_nwtrph(
+    void rhs_cellwise_nwtrph_ice(
       const opts_t<real_t> &opts,
       cont_t &dot_th_cont,
       cont_t &dot_rv_cont,
@@ -175,58 +142,109 @@ namespace libcloudphxx
       const cont_t &rc_cont,
       const cont_t &rr_cont,
       const cont_t &ria_cont,
-      const real_t &dt // time step in seconds
+      const real_t &dt
     )
     {
       rhs_cellwise<real_t, cont_t>(opts, dot_rc_cont, dot_rr_cont,  rc_cont, rr_cont);
-      if (opts.ice) {
-        rhs_cellwise_ice<real_t, cont_t>(opts, dot_rc_cont, dot_rr_cont, dot_rv_cont, dot_ria_cont, rc_cont, rr_cont, rv_cont, ria_cont, th_cont, p_cont, rhod_cont, dt);
-      }
 
       // rain evaporation treated as a force in Newthon-Raphson saturation adjustment
-      for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rr_cont, rhod_cont, p_cont, th_cont, rv_cont, rr_cont))
+      for (auto tup : zip(dot_th_cont, dot_rv_cont, dot_rc_cont, dot_rr_cont, dot_ria_cont, rhod_cont, p_cont, th_cont, rv_cont, rc_cont, rr_cont, ria_cont))
       {
         using namespace common;
 
         real_t
-          tmp = 0,
-          &dot_th = std::get<0>(tup),
-          &dot_rv = std::get<1>(tup),
-          &dot_rr = std::get<2>(tup);
+          rr_to_rv = 0,
+          rv_to_ria = 0,
+          rc_to_ria = 0,
+        &dot_th = std::get<0>(tup),
+        &dot_rv = std::get<1>(tup),
+        &dot_rc = std::get<2>(tup),
+        &dot_rr = std::get<3>(tup),
+        &dot_ria = std::get<4>(tup);
 
         const quantity<si::mass_density, real_t>
-          rhod = std::get<3>(tup) * si::kilograms / si::cubic_metres;
+          rhod = std::get<5>(tup) * si::kilograms / si::cubic_metres;
 
         const quantity<si::pressure, real_t>
-          p   = std::get<4>(tup) * si::pascals;
+          p   = std::get<6>(tup) * si::pascals;
 
         const quantity<si::temperature, real_t>
-          th = std::get<5>(tup) * si::kelvins;
+          th = std::get<7>(tup) * si::kelvins;
 
         const real_t
-          &rv = std::get<6>(tup),
-          &rr = std::get<7>(tup);
+          &rv     = std::get<8>(tup),
+          &rc     = std::get<9>(tup),
+          &rr     = std::get<10>(tup),
+          &ria     = std::get<11>(tup);
 
         quantity<si::temperature, real_t> T = th * theta_std::exner(p);
-        real_t r_vs = const_cp::r_vs(T, p);
-
-        tmp = (
-          formulae::evaporation_rate(
-            rv * si::dimensionless(),
-            r_vs * si::dimensionless(),
-            rr * si::dimensionless(),
-            rhod,
-            p
-          ) * si::seconds * dt
+        real_t rvs = const_cp::r_vs(T, p);
+        real_t rvsi = const_cp::r_vsi(T, p);
+
+        // rain evaporation
+        rr_to_rv += (
+        formulae::evaporation_rate(
+          rv * si::dimensionless(),
+          rvs * si::dimensionless(),
+          rr * si::dimensionless(),
+          rhod,
+          p
+        ) * si::seconds * dt
         );
 
         // limiting
-        tmp = std::min(rr, tmp) / dt;
+        rr_to_rv = std::min(rr, rr_to_rv) / dt;
 
-        dot_rv += tmp;
-        dot_rr -= tmp;
+        // ice A heterogeneous nucleation
+        if (opts.hetA)
+        {
+          rc_to_ria += (
+            formulae::het_A_nucleation(
+                    ria * si::dimensionless(),
+                    rc * si::dimensionless(),
+                    th,
+                    p,
+                    rhod,
+                    dt * si::seconds
+                  ) * si::seconds // to make it dimensionless
+          );
+        }
+
+        // ice A homogeneous nucleation 1
+        if (opts.homA1)
+        {
+          rv_to_ria += (
+            formulae::hom_A_nucleation_1(
+                    rv * si::dimensionless(),
+                    rvs * si::dimensionless(),
+                    rvsi * si::dimensionless(),
+                    th,
+                    p,
+                    dt * si::seconds
+                  ) * si::seconds // to make it dimensionless
+          );
+        }
+
+        // ice A homogeneous nucleation 2
+        if (opts.homA2)
+        {
+          rc_to_ria += (
+            formulae::hom_A_nucleation_2(
+                    rc * si::dimensionless(),
+                    th,
+                    p,
+                    dt * si::seconds
+                  ) * si::seconds // to make it dimensionless
+          );
+        }
 
-        dot_th -= const_cp::l_v(T) / (moist_air::c_pd<real_t>() * theta_std::exner(p)) * tmp / si::kelvins;
+        dot_rr -= rr_to_rv;
+        dot_rc -= rc_to_ria;
+        dot_rv += rr_to_rv - rv_to_ria;
+        dot_ria += rc_to_ria + rv_to_ria;
+        dot_th -= const_cp::l_v(T) / (moist_air::c_pd<real_t>() * theta_std::exner(p)) * rr_to_rv / si::kelvins; //heat of vaporisation
+        dot_th += const_cp::l_sublim(T) / (moist_air::c_pd<real_t>() * theta_std::exner(p)) * rv_to_ria / si::kelvins; //heat of sublimation
+        dot_th += const_cp::l_freez(T) / (moist_air::c_pd<real_t>() * theta_std::exner(p)) * rc_to_ria / si::kelvins; //heat of freezing
       }
     }
   }
diff --git a/include/libcloudph++/common/const_cp.hpp b/include/libcloudph++/common/const_cp.hpp
index 4fb387b0..71151795 100644
--- a/include/libcloudph++/common/const_cp.hpp
+++ b/include/libcloudph++/common/const_cp.hpp
@@ -23,6 +23,7 @@ namespace libcloudphxx
       libcloudphxx_const(si::temperature, T_tri, 273.16, si::kelvins) // temperature
       libcloudphxx_const(energy_over_mass, l_tri_evap, 2.5e6, si::joules / si::kilograms) // latent heat of evaporation
       libcloudphxx_const(energy_over_mass, l_tri_sublim, 2.83e6, si::joules / si::kilograms) // latent heat of sublimation
+      libcloudphxx_const(energy_over_mass, l_tri_freez, 3.34e5, si::joules / si::kilograms) // latent heat of freezing/melting
 
       // saturation vapour pressure with respect to liquid water
       // assuming constant c_p_v and c_p_w with constants taken at triple point
@@ -51,7 +52,6 @@ namespace libcloudphxx
       {
         return p_tri<real_t>() * exp(
           (l_tri_sublim<real_t>() + (c_pi<real_t>() - c_pv<real_t>()) * T_tri<real_t>()) / R_v<real_t>() * (real_t(1) / T_tri<real_t>() - real_t(1) / T)
-          - (c_pi<real_t>() - c_pv<real_t>()) / R_v<real_t>() * std::log(T / T_tri<real_t>())
         );
       }
 
@@ -83,6 +83,25 @@ namespace libcloudphxx
       ) {
         return l_tri_evap<real_t>() + (c_pv<real_t>() - c_pw<real_t>()) * (T - T_tri<real_t>());
       }
+
+      // latent heat of sublimation for constant c_p
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<divide_typeof_helper<si::energy, si::mass>::type , real_t> l_sublim(
+        const quantity<si::temperature, real_t> &T
+      ) {
+        return l_tri_sublim<real_t>() + (c_pv<real_t>() - c_pi<real_t>()) * (T - T_tri<real_t>());
+      }
+
+      // latent heat of freezing/melting for constant c_p
+      template <typename real_t>
+      BOOST_GPU_ENABLED
+      quantity<divide_typeof_helper<si::energy, si::mass>::type , real_t> l_freez(
+        const quantity<si::temperature, real_t> &T
+      ) {
+        return l_tri_freez<real_t>() + (c_pw<real_t>() - c_pi<real_t>()) * (T - T_tri<real_t>());
+      }
+
     };
   };
 };
diff --git a/tests/python/unit/api_blk_1m.py b/tests/python/unit/api_blk_1m.py
index 6e87a18d..1c36b53a 100644
--- a/tests/python/unit/api_blk_1m.py
+++ b/tests/python/unit/api_blk_1m.py
@@ -12,7 +12,6 @@
 print("conv =", opts.conv) 
 print("accr =", opts.accr) 
 print("sedi =", opts.sedi)
-print("ice =", opts.ice)
 print("homA1 =", opts.homA1)
 print("homA2 =", opts.homA2)
 print("hetA =", opts.hetA)
@@ -63,7 +62,6 @@
 
 dot_rc = arr_t([0.])
 dot_rr = arr_t([0.])
-dot_ria = arr_t([0.])
 blk_1m.rhs_cellwise(opts, dot_rc, dot_rr, rc, rr)
 assert dot_rc != 0 # some water should have coalesced
 assert dot_rr != 0
@@ -75,7 +73,7 @@
 dot_rc = arr_t([0.])
 dot_rr = arr_t([0.])
 
-blk_1m.rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, rhod, p, th, rv, rc, rr, ria, dt)
+blk_1m.rhs_cellwise_nwtrph(opts, dot_th, dot_rv, dot_rc, dot_rr, rhod, p, th, rv, rc, rr, dt)
 assert dot_rc != 0 # some water should have coalesced
 assert dot_rr != 0
 assert dot_th != 0 # some rain should have evaporated
@@ -92,6 +90,5 @@
 dot_rr = arr_t([0.])
 dot_rv = arr_t([0.])
 dot_ria = arr_t([0.])
-rv   = arr_t([0.01])
-blk_1m.rhs_cellwise_ice(opts, dot_rc, dot_rr, dot_rv, dot_ria, rc, rr, rv, ria, th, p, rhod, dt)
-assert dot_ria != 0
+blk_1m.rhs_cellwise_nwtrph_ice(opts, dot_th, dot_rv, dot_rc, dot_rr, dot_ria, rhod, p, th, rv, rc, rr, ria, dt)
+assert dot_ria != 0
\ No newline at end of file