diff --git a/py_wake/tests/test_sites/test_site.py b/py_wake/tests/test_sites/test_site.py
index 7ce1a9c433f7c14fb519cfbb7f16b84a1d78a4a8..77e358c6f51938886109c1f7f2c80a620ea2575c 100644
--- a/py_wake/tests/test_sites/test_site.py
+++ b/py_wake/tests/test_sites/test_site.py
@@ -16,6 +16,15 @@ k = [2.392578, 2.447266, 2.412109, 2.591797, 2.755859, 2.595703,
 ti = .1
 
 
+@pytest.fixture(autouse=True)
+def close_plots():
+    yield
+    try:
+        plt.close()
+    except Exception:
+        pass
+
+
 @pytest.fixture
 def site():
     return UniformWeibullSite(f, A, k, ti, shear=PowerShear(50, alpha=np.zeros_like(f) + .3))
diff --git a/py_wake/tests/test_sites/test_wasp_grid_site.py b/py_wake/tests/test_sites/test_wasp_grid_site.py
index 4930e4cc94723f7935050e9596899ff92b039659..1c2a75d5f9a8f57d9e26c8ee21b2f3b1cef576a6 100644
--- a/py_wake/tests/test_sites/test_wasp_grid_site.py
+++ b/py_wake/tests/test_sites/test_wasp_grid_site.py
@@ -16,6 +16,15 @@ from py_wake.tests.check_speed import timeit
 import matplotlib.pyplot as plt
 
 
+@pytest.fixture(autouse=True)
+def close_plots():
+    yield
+    try:
+        plt.close()
+    except Exception:
+        pass
+
+
 @pytest.fixture
 def site():
     return ParqueFicticioSite()
diff --git a/py_wake/tests/test_sites/test_xrsite.py b/py_wake/tests/test_sites/test_xrsite.py
index 83e76a57ad210c999ca8c9aa7ce26140f2d07c3d..7075e8677fe60051c1a71fecaef16560b8a38d70 100644
--- a/py_wake/tests/test_sites/test_xrsite.py
+++ b/py_wake/tests/test_sites/test_xrsite.py
@@ -27,6 +27,15 @@ k = [2.392578, 2.447266, 2.412109, 2.591797, 2.755859, 2.595703,
 ti = .1
 
 
+@pytest.fixture(autouse=True)
+def close_plots():
+    yield
+    try:
+        plt.close()
+    except Exception:
+        pass
+
+
 @pytest.fixture
 def uniform_site():
     ti = 0.1
diff --git a/py_wake/wind_farm_models/engineering_models.py b/py_wake/wind_farm_models/engineering_models.py
index df008bfffe3c92a1087b1ce3027b3ce12741a44b..b84b727ee35e9c2fefadf1e850dbeac9ed3d410d 100644
--- a/py_wake/wind_farm_models/engineering_models.py
+++ b/py_wake/wind_farm_models/engineering_models.py
@@ -272,6 +272,7 @@ class EngineeringWindFarmModel(WindFarmModel):
                     add_turb_ijk = self.turbulenceModel.calc_added_turbulence(dw_ijlk=dw_ijlk, **args)[:, :, 0]
 
             l_ = [l, 0][lw_j.WS_ilk.shape[1] == 1]
+            assert lw_j.WS_ilk.shape[1] == dh_ijl.shape[2] == hcw_ijlk.shape[2]
             if isinstance(self.superpositionModel, WeightedSum):
                 cw_ijk = np.hypot(dh_ijl[..., na], hcw_ijlk)[:, :, l_]
                 hcw_ijk, dh_ijk = hcw_ijlk[:, :, l_], dh_ijl[:, :, l_, na]