diff --git a/earthdiagnostics/ocean/heatcontentlayer.py b/earthdiagnostics/ocean/heatcontentlayer.py
index edf833260b8a800f53e816ef5191cb175a53b724..b5093e833f26d2e06ea89ed148c6a7cf34048cf9 100644
--- a/earthdiagnostics/ocean/heatcontentlayer.py
+++ b/earthdiagnostics/ocean/heatcontentlayer.py
@@ -47,7 +47,7 @@ class HeatContentLayer(Diagnostic):
     "Diagnostic alias for the configuration file"
 
     def __init__(self, data_manager, startdate, member, chunk, box, areas,
-                 weight, layers, basins, data_convention):
+                 weight, layers, basins, data_convention, min_level, max_level):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
         self.member = member
@@ -58,6 +58,8 @@ class HeatContentLayer(Diagnostic):
         self.layers = layers
         self.basins = basins
         self.data_convention = data_convention
+        self.min_level = min_level
+        self.max_level = max_level
 
     def __str__(self):
         return 'Heat content layer Startdate: {0} Member: {1} Chunk: {2} Box: {3}'.format(self.startdate, self.member,
@@ -106,6 +108,8 @@ class HeatContentLayer(Diagnostic):
         mask = mesh.get_landsea_mask()
         depth = mesh.get_depth(cell_point='W')
         weight = ohc.get_weights(layers, mask, e3t, depth)
+        max_level, min_level = cls._get_used_levels(weight)
+        weight[0] = weight[0][:, min_level:max_level,:, :]
 
         del mask, depth, e3t
 
@@ -113,7 +117,7 @@ class HeatContentLayer(Diagnostic):
             job_list.append(HeatContentLayer(
                 diags.data_manager, startdate, member, chunk, box,
                 areas, weight, layers, basins,
-                diags.config.data_convention
+                diags.config.data_convention, min_level, max_level
             ))
         return job_list
 
@@ -159,6 +163,18 @@ class HeatContentLayer(Diagnostic):
                                            self.startdate, self.member,
                                            self.chunk, box=self.box)
 
+    @classmethod
+    def _get_used_levels(cls, weight):
+        # Now we will reduce to the levels with any weight != 0 to avoid loading too much data on memory
+        levels = weight[0].shape[1]
+        min_level = 0
+        while min_level < levels and not weight[0][:, min_level, :].any():
+            min_level += 1
+        max_level = min_level
+        while max_level < (levels - 1) and weight[0][:, max_level + 1, :].any():
+            max_level += 1
+        return max_level, min_level
+
     def compute(self):
         """Run the diagnostic"""
         thetao_file = TempFile.get()
@@ -169,7 +185,7 @@ class HeatContentLayer(Diagnostic):
         handler = Utils.open_cdf(thetao_file)
         Utils.convert_units(handler.variables['thetao'], 'K')
         heatc_sl, heatc_sl1D = ohc.compute(self.layers, self.weight,
-                                           handler.variables['thetao'][:],
+                                           handler.variables['thetao'][:,self.min_level:self.max_level,:,:],
                                            self.areas)
         handler.sync()
         handler.renameVariable('thetao', 'heatc_sl')
@@ -213,3 +229,4 @@ class HeatContentLayer(Diagnostic):
         Utils.setminmax(results, 'heatc')
         self.heatc.set_local_file(results)
         self.heatcsum.set_local_file(results1D)
+