limit rhs to avoid negative values

src/solvers/blk_2m/slvr_blk_2m_common.hpp

@@ -76,12 +76,28 @@ class slvr_blk_2m_common : public std::conditional_t<ct_params_t::sgs_scheme ==
   void hook_mixed_rhs_ante_step()
   {
     negcheck(this->mem->advectee(ix::nc)(this->ijk), "nc at start of hook_mixed_rhs_ante_step");
-    //negtozero(this->mem->advectee(ix::nc)(this->ijk), "nc at start of hook_mixed_rhs_ante_step");
+  //  negtozero(this->mem->advectee(ix::nc)(this->ijk), "nc at start of hook_mixed_rhs_ante_step");
     parent_t::hook_mixed_rhs_ante_step();
+
+    nancheck(this->mem->advectee(ix::nc)(this->ijk), "nc after mixed_rhs_ante_step apply rhs");
+    nancheck(this->mem->advectee(ix::nr)(this->ijk), "nr after mixed_rhs_ante_step apply rhs");
+    nancheck(this->mem->advectee(ix::rc)(this->ijk), "rc after mixed_rhs_ante_step apply rhs");
+    nancheck(this->mem->advectee(ix::rr)(this->ijk), "rr after mixed_rhs_ante_step apply rhs");
+    negcheck2(this->mem->advectee(ix::nc)(this->ijk), this->rhs.at(ix::nc)(this->ijk), "nc after mixed_rhs_ante_step apply rhs (+ output of nc rhs)");
+//    negcheck2(this->mem->advectee(ix::nr)(this->ijk), this->rhs.at(ix::nr)(this->ijk), "nr after mixed_rhs_ante_step apply rhs (+ output of nr rhs)");
+    negcheck2(this->mem->advectee(ix::rc)(this->ijk), this->rhs.at(ix::rc)(this->ijk), "rc after mixed_rhs_ante_step apply rhs (+ output of rc rhs)");
+    //negcheck2(this->mem->advectee(ix::rr)(this->ijk), this->rhs.at(ix::rr)(this->ijk), "rr after mixed_rhs_ante_step apply rhs (+ output of rr rhs)");
+    negtozero(this->mem->advectee(ix::nr)(this->ijk), "nr after mixed_rhs_ante_step apply rhs");
+    negtozero(this->mem->advectee(ix::rr)(this->ijk), "rr after mixed_rhs_ante_step apply rhs");
   }
 
   void hook_mixed_rhs_post_step()
   {
+    negtozero(this->mem->advectee(ix::nc)(this->ijk), "nc at start of hook_mixed_rhs_ante_step");
+    negtozero(this->mem->advectee(ix::nr)(this->ijk), "nr at start of hook_mixed_rhs_ante_step");
+    negtozero(this->mem->advectee(ix::rr)(this->ijk), "rr at start of hook_mixed_rhs_ante_step");
+    //negcheck(this->mem->advectee(ix::nc)(this->ijk), "nc at start of mixed_rhs_post_step");
+
     parent_t::hook_mixed_rhs_post_step();
 
     nancheck(this->mem->advectee(ix::nc)(this->ijk), "nc after mixed_rhs_post_step apply rhs");
@@ -146,7 +162,7 @@ class slvr_blk_2m_common : public std::conditional_t<ct_params_t::sgs_scheme ==
     libmpdataxx::arrvec_t<typename parent_t::arr_t> &rhs,
     const typename parent_t::real_t &dt,
     const int &at
-  );
+  ) override;
 
   public:
 

src/solvers/blk_2m/update_rhs_blk_2m_common.hpp

@@ -52,6 +52,8 @@ void slvr_blk_2m_common<ct_params_t>::update_rhs(
       negcheck(nc, "nc after blk_2m rhs_cellwise call");
       nancheck(nc, "nr after blk_2m rhs_cellwise call");
       negcheck(nc, "nr after blk_2m rhs_cellwise call");
+
+//      negcheck(dot_nc, "dot_nc after blk_2m rhs_cellwise call");
   }
 
   // forcing
@@ -98,6 +100,27 @@ void slvr_blk_2m_common<ct_params_t>::update_rhs(
       break;
     }
   }
+
+  // assure that we do not remove more cloud/rain water than there is!
+  // this could happen due to blk_2m microphyscs combined with other forcings (e.g. subsidence)
+  auto
+    dot_rc = rhs.at(ix::rc)(this->ijk),
+    dot_rr = rhs.at(ix::rr)(this->ijk),
+    dot_nc = rhs.at(ix::nc)(this->ijk),
+    dot_nr = rhs.at(ix::nr)(this->ijk);
+
+  const auto
+    rc     = this->state(ix::rc)(this->ijk),
+    rr     = this->state(ix::rr)(this->ijk),
+    nc     = this->state(ix::nc)(this->ijk),
+    nr     = this->state(ix::nr)(this->ijk);
+
+  dot_rc = where(dot_rc * dt <= -rc, -rc / dt, dot_rc);
+  dot_rr = where(dot_rr * dt <= -rr, -rr / dt, dot_rr);
+  dot_nc = where(dot_nc * dt <= -nc, -nc / dt, dot_nc);
+  dot_nr = where(dot_nr * dt <= -rr, -rr / dt, dot_nr);
+
+
   this->mem->barrier();
   if(this->rank == 0)
   {