diff --git a/pydda/cost_functions/_cost_functions_numpy.py b/pydda/cost_functions/_cost_functions_numpy.py
index 95150117..9baf9c07 100644
--- a/pydda/cost_functions/_cost_functions_numpy.py
+++ b/pydda/cost_functions/_cost_functions_numpy.py
@@ -453,7 +453,7 @@ def calculate_mass_continuity(u, v, w, z, dx, dy, dz, coeff=1500.0, anel=1):
 
 
 def calculate_mass_continuity_gradient(
-    u, v, w, z, dx, dy, dz, coeff=1500.0, anel=1, upper_bc=True
+    u, v, w, z, dx, dy, dz, vrs=0, coeff=1500.0, anel=1, upper_bc=True, above=2.0
 ):
     """
     Calculates the gradient of mass continuity cost function. This is done by
@@ -500,10 +500,15 @@ def calculate_mass_continuity_gradient(
     grad_v = -np.gradient(div, dy, axis=1) * coeff
     grad_w = -np.gradient(div, dz, axis=0) * coeff
 
-    # Impermeability condition
-    grad_w[0, :, :] = 0
-    if upper_bc is True:
+    # Impermeability conditions
+    grad_w[0, :, :] = 0  # surface is impermeable
+    if upper_bc == 1:  # is True: # impermeable at the grid top
         grad_w[-1, :, :] = 0
+    if upper_bc == 2:  # impermeable at cloud top
+        N = np.sum(np.array(vrs) > -1000, axis=0)
+        z_mask = z > above * 1000
+        n_mask = N == 0
+        grad_w[z_mask & n_mask] = 0
     y = np.stack([grad_u, grad_v, grad_w], axis=0)
     return y.flatten()
 
diff --git a/pydda/cost_functions/cost_functions.py b/pydda/cost_functions/cost_functions.py
index 198e09b7..7787ffaf 100644
--- a/pydda/cost_functions/cost_functions.py
+++ b/pydda/cost_functions/cost_functions.py
@@ -543,9 +543,11 @@ def grad_J(winds, parameters):
                             parameters.dx,
                             parameters.dy,
                             parameters.dz,
+                            parameters.vrs,
                             parameters.Cm,
                             1,
                             parameters.upper_bc,
+                            above=parameters.above,
                         )
                     )
                 if parameters.Cx > 0 or parameters.Cy > 0 or parameters.Cz > 0:
diff --git a/pydda/retrieval/wind_retrieve.py b/pydda/retrieval/wind_retrieve.py
index ff22c37c..e757a8df 100644
--- a/pydda/retrieval/wind_retrieve.py
+++ b/pydda/retrieval/wind_retrieve.py
@@ -159,9 +159,12 @@ class DDParameters(object):
         Cartesian coordinates.
     roi: float
         The radius of influence of each point observation in m.
-    upper_bc: bool
-        True to enforce w=0 at top of domain (impermeability condition),
-        False to not enforce impermeability at top of domain
+    above: float
+        The altitude below which the cloud top impermeability would not apply (minimum impermeable height)
+    upper_bc: int
+        0 to not enforce impermeability at top of domain
+        1 to enforce w=0 at top of domain (impermeability condition),
+        2 to enforce w=0 at cloud top (or "above" if higher) (impermeability condition),
     """
 
     def __init__(self):
@@ -203,6 +206,7 @@ def __init__(self):
         self.upper_bc = True
         self.lower_bc = True
         self.roi = 1000.0
+        self.above = 2.0
         self.frz = 4500.0
         self.Nfeval = 0.0
         self.engine = "scipy"
@@ -254,10 +258,11 @@ def _get_dd_wind_field_scipy(
     leise_nstep=1,
     min_bca=30.0,
     max_bca=150.0,
-    upper_bc=True,
+    upper_bc=1,
     model_fields=None,
     output_cost_functions=True,
     roi=1000.0,
+    above=2.0,
     wind_tol=0.1,
     tolerance=1e-8,
     const_boundary_cond=False,