diff --git a/VERSION b/VERSION
index 2ea482c2b4b9d71cc91fa06f5a1c8c9fb373c1fb..bdb0fbcd421ce3d3d1524d61cf1a1c6e819dac3c 100644
--- a/VERSION
+++ b/VERSION
@@ -1,2 +1,2 @@
-3.0.0b50
+3.0.0b51
 
diff --git a/diags.conf b/diags.conf
index 360863d0e9d557fb0dd17e6cf3d658c046d51206..256eb0fc5a59d0d78b6791ab1404c932919078f7 100644
--- a/diags.conf
+++ b/diags.conf
@@ -15,7 +15,7 @@ CON_FILES = /esnas/autosubmit/con_files/
 # Diagnostics to run, space separated. You must provide for each one the name and the parameters (comma separated) or
 # an alias defined in the ALIAS section (see more below). If you are using the diagnostics just to CMORize, leave it
 # empty
-DIAGS = simdim,atmos,tas
+DIAGS = vgrad,thetao
 # DIAGS = OHC
 # Frequency of the data you want to use by default. Some diagnostics do not use this value: i.e. monmean always stores
 # its results at monthly frequency (obvious) and has a parameter to specify input's frequency.
@@ -85,7 +85,7 @@ OCEAN_TIMESTEP = 6
 # if 2, fc00
 # CHUNK_SIZE is the size of each data file, given in months
 # CHUNKS is the number of chunks. You can specify less chunks than present on the experiment
-EXPID = t01b
+EXPID = t01d
 STARTDATES = 19900101
 MEMBERS = 0
 MEMBER_DIGITS = 1
diff --git a/doc/source/codedoc/ocean.rst b/doc/source/codedoc/ocean.rst
index e8fcbe6e2bddf7e4909681c9df5d12053187cebf..cfa40e996173c8c5f98ae70c7ecca1da01900f13 100644
--- a/doc/source/codedoc/ocean.rst
+++ b/doc/source/codedoc/ocean.rst
@@ -103,6 +103,12 @@ earthdiagnostics.ocean.siasiesiv
    :show-inheritance:
    :members:
 
+earthdiagnostics.ocean.verticalgradient
+---------------------------------------
+.. automodule:: earthdiagnostics.ocean.verticalgradient
+   :show-inheritance:
+   :members:
+
 earthdiagnostics.ocean.verticalmean
 -----------------------------------
 .. automodule:: earthdiagnostics.ocean.verticalmean
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 4734cc4a185bd23c32a1243fcce014034af323f7..32893e405ba525e8255d0e990dc900e0f59d33fb 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -64,7 +64,7 @@ copyright = u'2016, BSC-CNS Earth Sciences Department'
 # The short X.Y version.
 version = '3.0b'
 # The full version, including alpha/beta/rc tags.
-release = '3.0.0b50'
+release = '3.0.0b51'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/doc/source/diagnostic_list.rst b/doc/source/diagnostic_list.rst
index f8f2d7efd93672c330ca9d8f5b016b93a9163f3b..88861569d0926cac799dc5ac5f2017029b7cb71e 100644
--- a/doc/source/diagnostic_list.rst
+++ b/doc/source/diagnostic_list.rst
@@ -594,6 +594,24 @@ Options:
 1. Basin = 'Global':
     Basin to restrict the computation to.
 
+vgrad
+~~~~~
+
+Calculates the gradient between two levels in a 3D ocean variable.
+See :class:`~earthdiagnostics.ocean.verticalgradient.VerticalGradient`
+
+Options:
+********
+
+1. Variable:
+    Variable to compute
+
+2. Upper level = 1:
+    Upper level. Will be used as the reference to compute the gradient
+
+3. Lower level = 2:
+    Lower level.
+
 verticalmean
 ~~~~~~~~~~~~
 
diff --git a/earthdiagnostics/EarthDiagnostics.pdf b/earthdiagnostics/EarthDiagnostics.pdf
index adda909f7f366adc6de17baf81e04b60908b6053..e789ce78e75869dcee7feb4ad83e61aebfc453f5 100644
Binary files a/earthdiagnostics/EarthDiagnostics.pdf and b/earthdiagnostics/EarthDiagnostics.pdf differ
diff --git a/earthdiagnostics/cmor_tables/default.csv b/earthdiagnostics/cmor_tables/default.csv
index 9224e2c95a3a517bf6aa1af7275bfc8e4f83de54..3a49dc44f2e50f79ec60f0d5b1044560822a3b07 100644
--- a/earthdiagnostics/cmor_tables/default.csv
+++ b/earthdiagnostics/cmor_tables/default.csv
@@ -343,4 +343,5 @@ zqla,hflso,surface_downward_latent_heat_flux,Surface Downward Latent Heat Flux,o
 zqsb,hfsso,surface_downward_sensible_heat_flux,Surface Downward Sensible Heat Flux,ocean,,W m-2,,,,
 zqlw,rlntds,surface_net_downward_longwave_flux,Surface Net Downward Longwave Radiation,ocean,,W m-2,,,,
 var78,tclw,total_column_liquid_water,Total column liquid water,atmos,,kg m-2,,,,
-var79,tciw,total_column_ice_water,Total column ice water,atmos,,kg m-2,,,,
\ No newline at end of file
+var79,tciw,total_column_ice_water,Total column ice water,atmos,,kg m-2,,,,
+rho,rhopoto,sea_water_potential_density,Sea Water Potential Density,ocean,,kg m-3,,,,
\ No newline at end of file
diff --git a/earthdiagnostics/cmor_tables/primavera b/earthdiagnostics/cmor_tables/primavera
index cbc0ea84e67caa713638e219d31a11dcc4b45039..ad4f256777265479b503bfa3e88a61ce05dd932f 160000
--- a/earthdiagnostics/cmor_tables/primavera
+++ b/earthdiagnostics/cmor_tables/primavera
@@ -1 +1 @@
-Subproject commit cbc0ea84e67caa713638e219d31a11dcc4b45039
+Subproject commit ad4f256777265479b503bfa3e88a61ce05dd932f
diff --git a/earthdiagnostics/cmorizer.py b/earthdiagnostics/cmorizer.py
index 76bbf50d93a08a0b4e1b8c393b404a8804a5c59d..b5590a6092b55a30ab9cc1bcf68fc7cf4f517d72 100644
--- a/earthdiagnostics/cmorizer.py
+++ b/earthdiagnostics/cmorizer.py
@@ -163,11 +163,9 @@ class Cmorizer(object):
         Utils.cdo.mergetime(input=gg_files, output=merged_gg)
         for filename in sh_files + gg_files:
             os.remove(filename)
-        Utils.nco.ncks(input=merged_sh, output=merged_gg, options='-A')
-        os.remove(merged_sh)
         tar_startdate = tarfile[0:-4].split('_')[5].split('-')
-        new_name = 'MMA_1m_{0[0]}_{0[1]}.nc'.format(tar_startdate)
-        shutil.move(merged_gg, os.path.join(self.cmor_scratch, new_name))
+        shutil.move(merged_gg, os.path.join(self.cmor_scratch, 'MMAGG_1m_{0[0]}_{0[1]}.nc'.format(tar_startdate)))
+        shutil.move(merged_sh, os.path.join(self.cmor_scratch, 'MMASH_1m_{0[0]}_{0[1]}.nc'.format(tar_startdate)))
 
     def cmorize_atmos(self):
         """
@@ -419,7 +417,7 @@ class Cmorizer(object):
 
     def get_date_str(self, file_path):
         file_parts = os.path.basename(file_path).split('_')
-        if file_parts[0] in (self.experiment.expid, 'MMA', 'MMO') or file_parts[0].startswith('ORCA'):
+        if file_parts[0] in (self.experiment.expid, 'MMA', 'MMASH', 'MMAGG', 'MMO') or file_parts[0].startswith('ORCA'):
             # Model output
             if file_parts[-1].endswith('.tar'):
                 file_parts = file_parts[-1][0:-4].split('-')
diff --git a/earthdiagnostics/cmormanager.py b/earthdiagnostics/cmormanager.py
index 13913fc8ffc9714a64f2445050f72d7e72198dc0..c34ffba8444642edfb89542e833489eec1044a32 100644
--- a/earthdiagnostics/cmormanager.py
+++ b/earthdiagnostics/cmormanager.py
@@ -474,8 +474,6 @@ class CMORManager(DataManager):
                 for filename in filenames:
                     if '_S{0}_'.format(startdate) in filename:
                         continue
-                    if self._get_member_str(member) in filename:
-                        continue
                     filepath = os.path.join(dirpath, filename)
                     good = filepath.replace('_{0}_output_'.format(self.experiment.model),
                                             '_{0}_{1}_S{2}_'.format(self.experiment.model,
@@ -487,7 +485,8 @@ class CMORManager(DataManager):
                                                           self.experiment.experiment_name))
 
                     Utils.move_file(filepath, good)
-                os.rmdir(dirpath)
+                if self.experiment.model != self.experiment.experiment_name:
+                    os.rmdir(dirpath)
             Log.debug('Done')
 
     def _remove_extra_output_folder(self):
diff --git a/earthdiagnostics/config.py b/earthdiagnostics/config.py
index ba1efb7ad244a6f55248092570c2735d92b677fd..7b339b61080c4b489b6533aa1be9cee06d697605 100644
--- a/earthdiagnostics/config.py
+++ b/earthdiagnostics/config.py
@@ -7,6 +7,7 @@ from bscearth.utils.config_parser import ConfigParser
 
 from earthdiagnostics.frequency import Frequency, Frequencies
 from earthdiagnostics.variable import VariableManager
+from modelingrealm import ModelingRealm
 
 
 class Config(object):
@@ -37,6 +38,7 @@ class Config(object):
         "Mask and meshes folder path"
         self.data_convention = parser.get_choice_option('DIAGNOSTICS', 'DATA_CONVENTION',
                                                         ('specs', 'primavera', 'cmip6'), 'specs', ignore_case=True)
+        VariableManager().load_variables(self.data_convention)
         self._diags = parser.get_option('DIAGNOSTICS', 'DIAGS')
         self.frequency = Frequency(parser.get_option('DIAGNOSTICS', 'FREQUENCY'))
         "Default data frequency to be used by the diagnostics"
@@ -106,10 +108,21 @@ class CMORConfig(object):
         self.source = parser.get_option('CMOR', 'SOURCE', 'to be filled')
 
         vars_string = parser.get_option('CMOR', 'VARIABLE_LIST', '')
+        var_manager = VariableManager()
         if vars_string:
             self._variable_list = list()
             for domain_var in vars_string.split(' '):
-                self._variable_list.append(domain_var.lower())
+                if domain_var.startswith('#'):
+                    break
+                splitted = domain_var.split(':')
+                cmor_var = var_manager.get_variable(splitted[1], silent=True)
+                if not cmor_var:
+                    continue
+                if ModelingRealm(splitted[0]) != cmor_var.domain:
+                    Log.warning('Domain {0} for variable {1} is not correct: is {2}', splitted[0], cmor_var.short_name,
+                                cmor_var.domain)
+                    return
+                self._variable_list.append('{0.domain}:{0.short_name}'.format(cmor_var))
         else:
             self._variable_list = None
 
@@ -129,7 +142,7 @@ class CMORConfig(object):
             return True
         if not var_cmor:
             return False
-        return '{0}:{1}'.format(var_cmor.domain, var_cmor.short_name).lower() in self._variable_list
+        return '{0}:{1}'.format(var_cmor.domain, var_cmor.short_name) in self._variable_list
 
     def any_required(self, variables):
         if self._variable_list is None:
diff --git a/earthdiagnostics/datamanager.py b/earthdiagnostics/datamanager.py
index bdc06a651488c3220425b9f191bb350bb4777e35..a946dd328cd4126fb3415306fcb6313b25d3a72f 100644
--- a/earthdiagnostics/datamanager.py
+++ b/earthdiagnostics/datamanager.py
@@ -28,7 +28,6 @@ class DataManager(object):
         self.experiment = config.experiment
         self._checked_vars = list()
         self.variable_list = VariableManager()
-        self.variable_list.load_variables(self.config.data_convention)
         UnitConversion.load_conversions()
         self.lock = threading.Lock()
 
@@ -217,7 +216,8 @@ class DataManager(object):
                 raise ValueError('Original file {0} does not exists'.format(filepath))
             if not os.path.isdir(os.path.dirname(link_path)):
                 Utils.create_folder_tree(os.path.dirname(link_path))
-            os.symlink(filepath, link_path)
+            relative_path = os.path.relpath(filepath, os.path.dirname(link_path))
+            os.symlink(relative_path, link_path)
         except:
             raise
         finally:
@@ -389,14 +389,14 @@ class NetCDFFile(object):
             handler.variables['lat'].short_name = 'lat'
             handler.variables['lat'].standard_name = 'latitude'
 
+    EQUIVALENT_UNITS = {'-': '1.0', 'fractional': '1.0', 'psu': 'psu'}
+
     def _fix_units(self, var_handler):
         if 'units' not in var_handler.ncattrs():
             return
-        if var_handler.units == '-':
-            var_handler.units = '1.0'
-        if var_handler.units == 'PSU':
-            var_handler.units = 'psu'
-        if var_handler.units == 'C' and self.cmor_var.units == 'K':
+        if var_handler.units.lower() in NetCDFFile.EQUIVALENT_UNITS:
+            var_handler.units = NetCDFFile.EQUIVALENT_UNITS[var_handler.units.lower()]
+        elif var_handler.units == 'C' or self.cmor_var.units == 'K':
             var_handler.units = 'deg_C'
         if self.cmor_var.units != var_handler.units:
             self._convert_units(var_handler)
diff --git a/earthdiagnostics/diagnostic.py b/earthdiagnostics/diagnostic.py
index 2a235977af263291d1826624c2e4d179d9a6e990..6d8036956f5e18758622b92102539c750054fb8e 100644
--- a/earthdiagnostics/diagnostic.py
+++ b/earthdiagnostics/diagnostic.py
@@ -219,6 +219,19 @@ class DiagnosticVariableOption(DiagnosticOption):
         return real_name.short_name
 
 
+class DiagnosticVariableListOption(DiagnosticOption):
+    def parse(self, option_value):
+        option_value = self.check_default(option_value)
+        var_names = []
+        for value in option_value.split('-'):
+            real_name = VariableManager().get_variable(value, False)
+            if real_name is None:
+                var_names.append(value)
+            else:
+                var_names.append(real_name.short_name)
+        return var_names
+
+
 class DiagnosticDomainOption(DiagnosticOption):
     def parse(self, option_value):
         return ModelingRealms.parse(self.check_default(option_value))
diff --git a/earthdiagnostics/earthdiags.py b/earthdiagnostics/earthdiags.py
index 5474f9b2215f119e25faf9c26b48627a8397004b..44fa9974496c4e38981eb8b112f5c02a905c2e6b 100755
--- a/earthdiagnostics/earthdiags.py
+++ b/earthdiagnostics/earthdiags.py
@@ -253,7 +253,7 @@ class EarthDiags(object):
         Diagnostic.register(RelinkAll)
         Diagnostic.register(Scale)
         Diagnostic.register(Attribute)
-        Diagnostic.register(SimplifyDimensions)
+        Diagnostic.register(SelectLevels)
 
     @staticmethod
     def _register_ocean_diagnostics():
@@ -277,6 +277,7 @@ class EarthDiags(object):
         Diagnostic.register(RegionMean)
         Diagnostic.register(Rotation)
         Diagnostic.register(Mxl)
+        Diagnostic.register(VerticalGradient)
 
     def clean(self):
         Log.info('Removing scratch folder...')
@@ -294,6 +295,7 @@ class EarthDiags(object):
                 report_path = os.path.join(self.config.scratch_dir,
                                            '{0}_{1}.report'.format(startdate,
                                                                    self.config.experiment.get_member_str(member)))
+                Utils.create_folder_tree(self.config.scratch_dir)
                 self.create_report(report_path, results)
 
         Log.result('Report finished')
@@ -312,28 +314,32 @@ class EarthDiags(object):
         return results
 
     def create_report(self, report_path, results):
-        current_table = None
-        current_priority = 0
-        results = sorted(results, key=lambda result: result[0].short_name)
-        results = sorted(results, key=lambda result: result[0].priority)
-        results = sorted(results, key=lambda result: result[1].name)
-
-        file_handler = open(report_path, 'w')
-
-        for var, table in results:
-            if current_table != table.name:
-                file_handler.write('\nTable {0}\n'.format(table.name))
+        realms = set([result[0].domain for result in results])
+        realms = sorted(realms)
+        for realm in realms:
+            file_handler = open('{0}.{1}'.format(report_path, realm), 'w')
+            realm_results = [result for result in results if result[0].domain == realm]
+
+            tables = set([result[1].name for result in realm_results])
+            tables = sorted(tables)
+            for table in tables:
+                table_results = [result for result in realm_results if result[1].name == table]
+
+                file_handler.write('\nTable {0}\n'.format(table))
                 file_handler.write('===================================\n')
-                current_table = table.name
-                current_priority = 0
 
-            if current_priority != var.priority:
-                file_handler.write('\nMissing variables with priority {0}:\n'.format(var.priority))
-                file_handler.write('--------------------------------------\n')
-                current_priority = var.priority
+                priorities = set([int(result[0].priority) for result in table_results])
+                priorities = sorted(priorities)
+                for priority in priorities:
+                    priority_results = [result for result in table_results if int(result[0].priority) == priority]
+                    priority_results = sorted(priority_results, key=lambda result: result[0].short_name)
+                    file_handler.write('\nMissing variables with priority {0}:\n'.format(priority))
+                    file_handler.write('--------------------------------------\n')
+
+                    for var, table in priority_results:
+                        file_handler.write('{0:12}: {1}\n'.format(var.short_name, var.standard_name))
 
-            file_handler.write('{0:12}: {1}\n'.format(var.short_name, var.standard_name))
-        file_handler.close()
+            file_handler.close()
 
     def _run_jobs(self, queue, numthread):
         def _run_job(current_job, retrials=1):
diff --git a/earthdiagnostics/frequency.py b/earthdiagnostics/frequency.py
index 35b19c2bfcb22c11a52de6accbcd8f645777b94a..1473fa5a2b8798923f405ba3ba5aaeb5086da53d 100644
--- a/earthdiagnostics/frequency.py
+++ b/earthdiagnostics/frequency.py
@@ -15,6 +15,7 @@ class Frequency(object):
                    '6': '6hr', '6h': '6hr', '6hr': '6hr', '6_hourly': '6hr', '6hourly': '6hr', '6 hourly': '6hr',
                    '3': '3hr', '3h': '3hr', '3hr': '3hr', '3_hourly': '3hr', '3hourly': '3hr', '3 hourly': '3hr',
                    '1': '1hr', 'hr': '1hr', '1h': '1hr', 'hourly': '1hr', '1hr': '1hr', '1 hourly': '1hr',
+                   '450mn': '450mn',
                    'subhr': 'subhr'}
 
     def __init__(self, freq):
diff --git a/earthdiagnostics/general/__init__.py b/earthdiagnostics/general/__init__.py
index 96be392dee4fb6951f742899d15767baa0b5deab..2920230a23576364f5a0b6d12559a6c8f6966e56 100644
--- a/earthdiagnostics/general/__init__.py
+++ b/earthdiagnostics/general/__init__.py
@@ -8,3 +8,4 @@ from earthdiagnostics.general.scale import Scale
 from earthdiagnostics.general.attribute import Attribute
 from earthdiagnostics.general.relinkall import RelinkAll
 from earthdiagnostics.general.simplify_dimensions import SimplifyDimensions
+from earthdiagnostics.general.select_levels import SelectLevels
diff --git a/earthdiagnostics/general/select_levels.py b/earthdiagnostics/general/select_levels.py
new file mode 100644
index 0000000000000000000000000000000000000000..a1fe6d0b01fe8d5193c0d7a0f0ff4a6f9168dc9a
--- /dev/null
+++ b/earthdiagnostics/general/select_levels.py
@@ -0,0 +1,111 @@
+# coding=utf-8
+from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, DiagnosticDomainOption, \
+    DiagnosticVariableListOption, DiagnosticIntOption
+from earthdiagnostics.modelingrealm import ModelingRealm
+from earthdiagnostics.utils import Utils, TempFile
+from earthdiagnostics.box import Box
+
+
+class SelectLevels(Diagnostic):
+    """
+    Convert i j files to lon lat when there is no interpolation required,
+    i.e. lon is constant over i and lat is constat over j
+
+    :original author: Javier Vegas-Regidor<javier.vegas@bsc.es>
+
+    :created: April 2017
+
+    :param data_manager: data management object
+    :type data_manager: DataManager
+    :param startdate: startdate
+    :type startdate: str
+    :param member: member number
+    :type member: int
+    :param chunk: chunk's number
+    :type chunk: int
+    :param variable: variable's name
+    :type variable: str
+    :param domain: variable's domain
+    :type domain: ModelingRealm
+    """
+
+    alias = 'selev'
+    "Diagnostic alias for the configuration file"
+
+    def __init__(self, data_manager, startdate, member, chunk, domain, variable, grid, first_level, last_level):
+        Diagnostic.__init__(self, data_manager)
+        self.startdate = startdate
+        self.member = member
+        self.chunk = chunk
+        self.variable = variable
+        self.domain = domain
+        self.grid = grid
+        self.box = Box(False)
+        self.box.min_depth = first_level
+        self.box.max_depth = last_level
+
+    def __str__(self):
+        return 'Select levels Startdate: {0} Member: {1} Chunk: {2} ' \
+               'Variable: {3}:{4} Levels: {6}-{7} Grid: {5}'.format(self.startdate, self.member, self.chunk, self.domain, self.variable,
+                                                    self.grid, self.box.min_depth, self.box.max_depth)
+
+    def __eq__(self, other):
+        return self.startdate == other.startdate and self.member == other.member and self.chunk == other.chunk and \
+            self.domain == other.domain and self.variable == other.variable and self.grid == self.grid
+
+    @classmethod
+    def generate_jobs(cls, diags, options):
+        """
+            Creates a job for each chunk to compute the diagnostic
+
+            :param diags: Diagnostics manager class
+            :type diags: Diags
+            :param options: domain,variables,grid
+            :type options: list[str]
+            :return:
+            """
+        options_available = (DiagnosticDomainOption('domain'),
+                             DiagnosticVariableListOption('variables'),
+                             DiagnosticIntOption('first_level'),
+                             DiagnosticIntOption('last_level'),
+                             DiagnosticOption('grid', ''))
+        options = cls.process_options(options, options_available)
+        job_list = list()
+        variables = options['variables']
+        for var in variables:
+            for startdate, member, chunk in diags.config.experiment.get_chunk_list():
+
+                job_list.append(SelectLevels(diags.data_manager, startdate, member, chunk,
+                                             options['domain'], var, options['grid'],
+                                             options['first_level'],options['last_level']))
+        return job_list
+
+    def compute(self):
+        """
+        Runs the diagnostic
+        """
+        variable_file = self.data_manager.get_file(self.domain, self.variable, self.startdate, self.member, self.chunk,
+                                                   grid=self.grid)
+
+        Utils.nco.ncks(input=variable_file, output=variable_file, options='-O -d lev,{0.min_depth},{0.max_depth}'.format(self.box))
+        self.send_file(variable_file, self.domain, self.variable, self.startdate, self.member, self.chunk,
+                       grid=self.grid)
+
+    def _create_var(self, var_name, var_values, source, destiny):
+        old_var = source.variables[var_name]
+        new_var = destiny.createVariable(var_name, old_var.dtype, dimensions=(var_name, ))
+        new_var[:] = var_values
+        Utils.copy_attributes(new_var, old_var)
+
+        vertices_name = '{0}_vertices'.format(var_name)
+
+        if vertices_name in source.variables:
+            var_vertices = source.variables[vertices_name]
+            if var_name == 'lon':
+                vertices_values = var_vertices[0:1, ...]
+            else:
+                vertices_values = var_vertices[:, 0:1, :]
+            new_lat_vertices = destiny.createVariable(vertices_name, var_vertices.dtype, dimensions=(var_name, 'vertices'))
+            new_lat_vertices[:] = vertices_values
+            Utils.copy_attributes(new_lat_vertices, var_vertices)
+
diff --git a/earthdiagnostics/general/simplify_dimensions.py b/earthdiagnostics/general/simplify_dimensions.py
index aeb76e4330403f916f2d29eff1e54bca1b35940b..63444c1da8c24e7452b5c49534beb4667abfb182 100644
--- a/earthdiagnostics/general/simplify_dimensions.py
+++ b/earthdiagnostics/general/simplify_dimensions.py
@@ -1,5 +1,6 @@
 # coding=utf-8
-from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, DiagnosticDomainOption, DiagnosticVariableOption
+from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, DiagnosticDomainOption, \
+    DiagnosticVariableListOption
 from earthdiagnostics.modelingrealm import ModelingRealm
 from earthdiagnostics.utils import Utils, TempFile
 import numpy as np
@@ -56,18 +57,21 @@ class SimplifyDimensions(Diagnostic):
 
             :param diags: Diagnostics manager class
             :type diags: Diags
-            :param options: variable, domain, grid
+            :param options: domain,variables,grid
             :type options: list[str]
             :return:
             """
         options_available = (DiagnosticDomainOption('domain'),
-                             DiagnosticVariableOption('variable'),
+                             DiagnosticVariableListOption('variables'),
                              DiagnosticOption('grid', ''))
         options = cls.process_options(options, options_available)
         job_list = list()
-        for startdate, member, chunk in diags.config.experiment.get_chunk_list():
-            job_list.append(SimplifyDimensions(diags.data_manager, startdate, member, chunk,
-                                               options['domain'], options['variable'], options['grid']))
+        variables = options['variables']
+        for var in variables:
+            for startdate, member, chunk in diags.config.experiment.get_chunk_list():
+
+                job_list.append(SimplifyDimensions(diags.data_manager, startdate, member, chunk,
+                                                   options['domain'], var, options['grid']))
         return job_list
 
     def compute(self):
diff --git a/earthdiagnostics/modelingrealm.py b/earthdiagnostics/modelingrealm.py
index acc5646c9c4b8b999dd200725f83e1e4bbf7ac90..caecd42d5605eea47267118f8308c93e1fcdd2dd 100644
--- a/earthdiagnostics/modelingrealm.py
+++ b/earthdiagnostics/modelingrealm.py
@@ -22,6 +22,9 @@ class ModelingRealm(object):
     def __eq__(self, other):
         return other.__class__ == ModelingRealm and self.name == other.name
 
+    def __ne__(self, other):
+        return not (self == other)
+
     def __str__(self):
         return self.name
 
diff --git a/earthdiagnostics/ocean/__init__.py b/earthdiagnostics/ocean/__init__.py
index 30db8957001262e3aa3c5ae0a7779eb9b0d70984..7512b3610c4c7db363a4fdf5e5ca1dd8fcc9b4bd 100644
--- a/earthdiagnostics/ocean/__init__.py
+++ b/earthdiagnostics/ocean/__init__.py
@@ -22,3 +22,4 @@ from earthdiagnostics.ocean.mixedlayerheatcontent import MixedLayerHeatContent
 from earthdiagnostics.ocean.regionmean import RegionMean
 from earthdiagnostics.ocean.rotation import Rotation
 from earthdiagnostics.ocean.mxl import Mxl
+from earthdiagnostics.ocean.verticalgradient import VerticalGradient
diff --git a/earthdiagnostics/ocean/heatcontent.py b/earthdiagnostics/ocean/heatcontent.py
index c67555828cf9f87e88258080a5e7104a5ac1c466..39729a9a62806d20bf98fc18dacfbf1ddaf44ac6 100644
--- a/earthdiagnostics/ocean/heatcontent.py
+++ b/earthdiagnostics/ocean/heatcontent.py
@@ -7,6 +7,7 @@ from earthdiagnostics.utils import Utils, TempFile
 from earthdiagnostics.diagnostic import Diagnostic, DiagnosticBasinOption, DiagnosticIntOption
 from earthdiagnostics.box import Box
 from earthdiagnostics.modelingrealm import ModelingRealms
+import numpy as np
 
 
 class HeatContent(Diagnostic):
@@ -37,7 +38,7 @@ class HeatContent(Diagnostic):
     alias = 'ohc'
     "Diagnostic alias for the configuration file"
 
-    def __init__(self, data_manager, startdate, member, chunk, basin, mixed_layer, box):
+    def __init__(self, data_manager, startdate, member, chunk, basin, mixed_layer, box, min_level, max_level):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
         self.member = member
@@ -45,6 +46,8 @@ class HeatContent(Diagnostic):
         self.basin = basin
         self.mxloption = mixed_layer
         self.box = box
+        self.min_level = min_level
+        self.max_level = max_level
         self.required_vars = ['so', 'mlotst']
         self.generated_vars = ['scvertsum']
 
@@ -74,13 +77,41 @@ class HeatContent(Diagnostic):
                              DiagnosticIntOption('min_depth'),
                              DiagnosticIntOption('max_depth'))
         options = cls.process_options(options, options_available)
-        box = Box(False)
+        box = Box(True)
         box.min_depth = options['min_depth']
         box.max_depth = options['max_depth']
+        min_level = 0
+        max_level = 0
+
+        if box.min_depth > 0 or box.max_depth > 0:
+
+            handler = Utils.openCdf('mesh_zgr.nc')
+
+            if 'gdepw_1d' in handler.variables:
+                depth = handler.variables['gdepw_1d'][0, :]
+            else:
+                raise Exception('gdepw 1D variable can not be found')
+            handler.close()
+
+            def find_nearest(array, value):
+                idx = (np.abs(array - value)).argmin()
+                return idx, round(array[idx])
+
+            if box.min_depth > 0:
+                min_level, box.min_depth = find_nearest(depth, box.min_depth)
+            else:
+                min_level = 1
+
+            if box.max_depth > 0:
+                max_level, box.max_depth = find_nearest(depth, box.max_depth)
+            else:
+                max_level = len(depth)
+                box.max_depth = round(depth[-1])
+
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job_list.append(HeatContent(diags.data_manager, startdate, member, chunk,
-                                        options['basin'], options['mixed_layer'], box))
+                                        options['basin'], options['mixed_layer'], box, min_level, max_level))
         return job_list
 
     def compute(self):
@@ -102,8 +133,8 @@ class HeatContent(Diagnostic):
             para.append(0)
             para.append(0)
             para.append(0)
-            para.append(self.box.min_depth)
-            para.append(self.box.max_depth)
+            para.append(self.min_level)
+            para.append(self.max_level)
         if self.mxloption != 0:
             para.append('-mxloption')
             para.append(str(self.mxloption))
diff --git a/earthdiagnostics/ocean/verticalgradient.py b/earthdiagnostics/ocean/verticalgradient.py
new file mode 100644
index 0000000000000000000000000000000000000000..1893aac892f4b0f04ce0e14e7af4793d881b98f2
--- /dev/null
+++ b/earthdiagnostics/ocean/verticalgradient.py
@@ -0,0 +1,114 @@
+# coding=utf-8
+from earthdiagnostics import cdftools
+from earthdiagnostics.box import Box
+from earthdiagnostics.diagnostic import Diagnostic, DiagnosticIntOption, DiagnosticVariableOption
+from earthdiagnostics.utils import Utils, TempFile
+from earthdiagnostics.modelingrealm import ModelingRealms
+
+
+class VerticalGradient(Diagnostic):
+    """
+    Chooses vertical level in ocean, or vertically averages between
+    2 or more  ocean levels
+
+    :original author: Virginie Guemas <virginie.guemas@bsc.es>
+    :contributor: Eleftheria Exarchou <eleftheria.exarchou@bsc.es>
+    :contributor: Javier Vegas-Regidor <javier.vegas@bsc.es>
+
+    :created: February 2012
+    :last modified: June 2016
+
+    :param data_manager: data management object
+    :type data_manager: DataManager
+    :param startdate: startdate
+    :type startdate: str
+    :param member: member number
+    :type member: int
+    :param chunk: chunk's number
+    :type chunk: int
+    :param variable: variable to average
+    :type variable: str
+    :param box: box used to restrict the vertical mean
+    :type box: Box
+    """
+
+    alias = 'vgrad'
+    "Diagnostic alias for the configuration file"
+
+    def __init__(self, data_manager, startdate, member, chunk, variable, box):
+        Diagnostic.__init__(self, data_manager)
+        self.startdate = startdate
+        self.member = member
+        self.chunk = chunk
+        self.variable = variable
+        self.box = box
+        self.required_vars = [variable]
+        self.generated_vars = [variable + 'vmean']
+
+    def __eq__(self, other):
+        return self.startdate == other.startdate and self.member == other.member and self.chunk == other.chunk and \
+            self.box == other.box and self.variable == other.variable
+
+    def __str__(self):
+        return 'Vertical gradient Startdate: {0} Member: {1} Chunk: {2} Variable: {3} ' \
+               'Box: {4}'.format(self.startdate, self.member, self.chunk, self.variable, self.box)
+
+    @classmethod
+    def generate_jobs(cls, diags, options):
+        """
+        Creates a job for each chunk to compute the diagnostic
+
+        :param diags: Diagnostics manager class
+        :type diags: Diags
+        :param options: variable, minimum depth (level), maximum depth (level)
+        :type options: list[str]
+        :return:
+        """
+        options_available = (DiagnosticVariableOption('variable'),
+                             DiagnosticIntOption('upper_level', 1),
+                             DiagnosticIntOption('low_level', 2))
+        options = cls.process_options(options, options_available)
+
+        box = Box(False)
+        if options['upper_level'] >= 0:
+            box.min_depth = options['upper_level']
+        if options['low_level'] >= 0:
+            box.max_depth = options['low_level']
+
+        job_list = list()
+        for startdate, member, chunk in diags.config.experiment.get_chunk_list():
+            job_list.append(VerticalGradient(diags.data_manager, startdate, member, chunk,
+                                             options['variable'], box))
+        return job_list
+
+    def compute(self):
+        """
+        Runs the diagnostic
+        """
+        variable_file = self.data_manager.get_file(ModelingRealms.ocean, self.variable, self.startdate, self.member,
+                                                   self.chunk)
+
+        handler = Utils.openCdf(variable_file)
+        if 'lev' not in handler.dimensions:
+            raise Exception('Variable {0} does not have a level dimension')
+        var_handler = handler.variables[self.variable]
+        upper_level = var_handler[:, self.box.min_depth-1, ...]
+        lower_level = var_handler[:, self.box.max_depth-1, ...]
+        gradient = upper_level - lower_level
+
+        temp = TempFile.get()
+        new_file = Utils.openCdf(temp, 'w')
+        for var in handler.variables.keys():
+            if var in (self.variable, 'lev', 'lev_bnds'):
+                continue
+            Utils.copy_variable(handler, new_file, var, add_dimensions=True)
+        new_var = new_file.createVariable(self.variable + 'vgrad',  var_handler.dtype,
+                                          dimensions=('time', 'j', 'i'), zlib=True)
+        Utils.copy_attributes(new_var, var_handler)
+        new_var[...] = gradient[...]
+        new_var.long_name += ' Vertical gradient'
+        new_var.standard_name += '_vertical_gradient'
+
+        self.send_file(temp, ModelingRealms.ocean, self.variable + 'vgrad', self.startdate, self.member, self.chunk,
+                       box=self.box)
+
diff --git a/earthdiagnostics/variable.py b/earthdiagnostics/variable.py
index 0ad73fba648335ef21902226dbdfb5cad7779349..d0bce08daeebed94870e7605a45bf0f268a24fbe 100644
--- a/earthdiagnostics/variable.py
+++ b/earthdiagnostics/variable.py
@@ -1,5 +1,6 @@
 # coding=utf-8
 import csv
+import glob
 import json
 
 import openpyxl
@@ -89,16 +90,16 @@ class VariableManager(object):
 
     def _load_variable_list(self, table_name):
 
+        xlsx_path = self._get_xlsx_path(table_name)
+        if xlsx_path:
+            self._load_xlsx(xlsx_path)
+            return
+
         json_folder = self._get_json_folder(table_name)
         if os.path.isdir(json_folder):
             self._load_json(json_folder)
             return
 
-        xlsx_path = self._get_xlsx_path(table_name)
-        if os.path.isfile(xlsx_path):
-            self._load_xlsx(table_name)
-            return
-
         csv_path = self._get_csv_path(table_name)
         if os.path.isfile(csv_path):
             self._load_file(table_name)
@@ -223,49 +224,58 @@ class VariableManager(object):
 
     def _get_xlsx_path(self, table_name):
         xlsx_table_path = os.path.join(self._cmor_tables_folder, '{0}.xlsx'.format(table_name))
-        return xlsx_table_path
 
-    def _load_xlsx(self, table_name):
-        xlsx_table_path = os.path.join(self._cmor_tables_folder, '{0}.xlsx'.format(table_name))
+        if os.path.isfile(xlsx_table_path):
+            return xlsx_table_path
+        xlsx_table_path = os.path.join(self._cmor_tables_folder, table_name, 'etc', '*.xlsx')
+        xlsx_table_path = glob.glob(xlsx_table_path)
+        if len(xlsx_table_path) == 1:
+            return xlsx_table_path[0]
+        return None
+
+    def _load_xlsx(self, xlsx_table_path):
         excel = openpyxl.load_workbook(xlsx_table_path, True)
 
         table_data = {}
         data_sheet = excel.worksheets[0]
         for row in data_sheet.rows:
             if row[1].value in excel.sheetnames:
-                table_data[row[1].value] = (Frequency(row[2].value), row[4].value[4:-1])
+                table_data[row[1].value] = (Frequency(row[2].value), 'Date missing')
 
         for sheet_name in excel.sheetnames:
-            sheet = excel.get_sheet_by_name(sheet_name)
-            if sheet['A1'].value != 'Priority':
-                continue
-            table_frequency, table_date = table_data[sheet.title]
-            table = CMORTable(sheet.title, table_frequency, table_date)
-            for row in sheet.rows:
-                if row[0].value == 'Priority' or not row[5].value:
+            try:
+                sheet = excel.get_sheet_by_name(sheet_name)
+                if sheet['A1'].value != 'Priority':
                     continue
+                table_frequency, table_date = table_data[sheet.title]
+                table = CMORTable(sheet.title, table_frequency, table_date)
+                for row in sheet.rows:
+                    if row[0].value == 'Priority' or not row[5].value:
+                        continue
 
-                if row[5].value.lower() in self._dict_variables:
-                    self._dict_variables[row[5].value.lower()].tables.append(table)
-                    continue
+                    if row[5].value.lower() in self._dict_variables:
+                        self._dict_variables[row[5].value.lower()].tables.append(table)
+                        continue
 
-                var = Variable()
-                var.priority = row[0].value
-                var.short_name = row[5].value
-                var.standard_name = row[6].value
-                var.long_name = row[1].value
+                    var = Variable()
+                    var.priority = row[0].value
+                    var.short_name = row[5].value
+                    var.standard_name = row[6].value
+                    var.long_name = row[1].value
 
-                self._process_modelling_realm(var, row[12].value)
+                    var.domain = self._process_modelling_realm(var, row[12].value)
 
-                var.units = row[2].value
-                var.tables.append(table)
-                self._dict_variables[var.short_name.lower()] = var
+                    var.units = row[2].value
+                    var.tables.append(table)
+                    self._dict_variables[var.short_name.lower()] = var
+            except Exception as ex:
+                Log.error('Table {0} can not be loaded: {1}', sheet_name, ex)
 
     def _process_modelling_realm(self, var, value):
         if value is None:
             value = ''
         modelling_realm = value.split(' ')
-        var.get_modelling_realm(modelling_realm)
+        return var.get_modelling_realm(modelling_realm)
 
     def _load_missing_defaults(self):
         self._load_file('default', True)
@@ -280,6 +290,9 @@ class Variable(object):
     def __str__(self):
         return '{0} ({1})'.format(self.standard_name, self.short_name)
 
+    def __repr__(self):
+        return '{0} ({1})'.format(self.standard_name, self.short_name)
+
     def __init__(self):
         self.short_name = None
         self.standard_name = None
@@ -309,6 +322,12 @@ class Variable(object):
         self.valid_min = json_var['valid_min'].strip()
         self.valid_max = json_var['valid_max'].strip()
         self.units = json_var['units'].strip()
+        if 'priority' in json_var:
+            self.priority = int(json_var['priority'].strip())
+        elif 'primavera_priority' in json_var:
+            self.priority = int(json_var['primavera_priority'].strip())
+        else:
+            self.priority = 1
 
     def get_modelling_realm(self, domains):
         if len(domains) > 1:
@@ -350,8 +369,10 @@ class Variable(object):
         for table in self.tables:
             if table.frequency == frequency:
                 return table
-        table_name = self.domain.get_table_name(frequency, data_convention)
-        return CMORTable(table_name, frequency, 'December 2013')
+        if self.domain:
+            table_name = self.domain.get_table_name(frequency, data_convention)
+            return CMORTable(table_name, frequency, 'December 2013')
+        return self.tables[0]
 
     def _select_most_specific(self, parsed):
         parsed = set(parsed)
@@ -399,4 +420,11 @@ class CMORTable(object):
     def __str__(self):
         return self.name
 
+    def __repr__(self):
+        return '{0.name} ({0.frequency}, {0.date}'.format(self)
+
+    def __lt__(self, other):
+        return self.name < other.name
+
+
 
diff --git a/earthdiagnostics/variable_alias/default.csv b/earthdiagnostics/variable_alias/default.csv
index 4a46d64d9832fbb0ac17699fad9d5d4a7054d1a5..8b967a4a575f19f1bde49a83f803a411ea4c8dc3 100644
--- a/earthdiagnostics/variable_alias/default.csv
+++ b/earthdiagnostics/variable_alias/default.csv
@@ -296,4 +296,5 @@ difvho,difvho,,
 vovematr,wmo,,
 qtr_ice,qtr,,
 var78,tclw,,
-var79,tciw,,
\ No newline at end of file
+var79,tciw,,
+rho,rhopoto,,
\ No newline at end of file
diff --git a/model_diags.conf b/model_diags.conf
index b55c2bc026459ee0da1c4f1b60e85e36ff0b043e..d479e349cf0f790bc4db223e985106b5aecdad6b 100644
--- a/model_diags.conf
+++ b/model_diags.conf
@@ -166,9 +166,9 @@ LMSALC = vertmeanmeters,so,300,5400
 USALC = vertmeanmeters,so,0,300
 OHC = ohc,glob,0,1,10
 XOHC = ohc,glob,1,0,0
-LOHC = ohc,glob,0,23,46
-MOHC = ohc,glob,0,18,22
-UOHC = ohc,glob,0,1,17
+LOHC = ohc,glob,0,700,2000
+MOHC = ohc,glob,0,300,700
+UOHC = ohc,glob,0,0,300
 OHC_SPECIFIED_LAYER = ohclayer,0,300 ohclayer,300,800
 3DTEMP = interp,thetao
 3DSAL = interp,so