diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index e7a4660d91d09fe1e621d0b8bec1266274149b3e..f3e8d8c5d078774ad1fd0a7eef79ba296be9dfc9 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -40,6 +40,7 @@ test_python3:
 report_codacy:
   stage: report
   script:
+    - '[ -z "$CODACY_PROJECT_TOKEN" ] && exit 0'
     - source activate earthdiagnostics3
     - pip install codacy-coverage --upgrade
     - python-codacy-coverage -r test/report/python3/coverage.xml
diff --git a/earthdiagnostics/cmor_tables/default.csv b/earthdiagnostics/cmor_tables/default.csv
index 7133c04f170f0df4e627a9082561e595ea417629..2fe501149d0fdffd78f5d53d76e569075bef4d34 100644
--- a/earthdiagnostics/cmor_tables/default.csv
+++ b/earthdiagnostics/cmor_tables/default.csv
@@ -1,5 +1,5 @@
 Variable,Shortname,Name,Long name,Domain,Basin,Units,Valid min,Valid max,Grid,Tables
-hfxout,hfxout,,Ocean surface heat flux,ocean,,W m-2,,,,,
+hfxout,hfxout,,Ocean surface heat flux,ocean,,W m-2,,,,
 iiceages:siage:iice_otd,ageice,age_of_sea_ice,Age of sea ice,seaIce,,,,,,
 al,al,surface_albedo,Albedo,atmos,,,,,,
 bgfrcsal,bgfrcsal,change_over_time_in_heat_content_from_forcing,Change over time in salt content from forcing,ocean,,,,,,
@@ -347,4 +347,7 @@ zqlw,rlntds,surface_net_downward_longwave_flux,Surface Net Downward Longwave Rad
 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,,,,
 rho,rhopoto,sea_water_potential_density,Sea Water Potential Density,ocean,,kg m-3,,,,
-hc300,hc300,,Heat Content to 300m depth,ocean,J m-2,,,,
\ No newline at end of file
+hc300,hc300,,Heat Content to 300m depth,ocean,,J m-2,,,,
+hfcorr,hfcorr,heat_flux_into_sea_water_due_to_newtonian_relaxation,heat_flux_correction,ocean,,W m-2,,,,
+wfcorr,wfcorr,water_flux_out_of_sea_water_due_to_newtonian_relaxation,water_flux_correction,ocean,,kg m-2 s-1,,,,
+hfxin,hfxin,,total heat fluxes at the ice/ocean surface,ocean,,W m-2,,,,
diff --git a/earthdiagnostics/config.py b/earthdiagnostics/config.py
index 6bdf281a20c0b905aabbef2c9f9ea970b8acd5a4..838bb223d5badda52511a04f6488986fc1dbc817 100644
--- a/earthdiagnostics/config.py
+++ b/earthdiagnostics/config.py
@@ -250,6 +250,7 @@ class CMORConfig(object):
         self.activity = parser.get_option('CMOR', 'ACTIVITY', 'CMIP')
         self.min_cmorized_vars = parser.get_int_option('CMOR', 'MIN_CMORIZED_VARS', 10)
         self.append_startdate = parser.get_bool_option('CMOR', 'APPEND_STARTDATE', False)
+        self.append_startdate_year_only = parser.get_bool_option('CMOR', 'APPEND_STARTDATE_YEAR_ONLY', False)
 
         vars_string = parser.get_option('CMOR', 'VARIABLE_LIST', '')
         self.var_manager = var_manager
diff --git a/earthdiagnostics/data_convention.py b/earthdiagnostics/data_convention.py
index 9e855dd16d880f645187ec661272fe9f4ec1695f..1c54d35cca874dc2430a1525eb285c4d09e5fa9b 100644
--- a/earthdiagnostics/data_convention.py
+++ b/earthdiagnostics/data_convention.py
@@ -585,6 +585,8 @@ class Cmor3Convention(DataConvention):
             else:
                 grid = self.config.cmor.default_atmos_grid
         if self.config.cmor.append_startdate:
+            if self.config.cmor.append_startdate_year_only:
+                startdate = startdate[0:4]
             subexp_id = 's{}-'.format(startdate)
         else:
             subexp_id = ''
diff --git a/earthdiagnostics/datafile.py b/earthdiagnostics/datafile.py
index a3d0ad22ea89153d8ab96fe6d5d8d80e8de44464..7c903e1e057a4f538b0c060a450ff55ea4eca40b 100644
--- a/earthdiagnostics/datafile.py
+++ b/earthdiagnostics/datafile.py
@@ -648,11 +648,11 @@ class NetCDFFile(DataFile):
             Log.debug('Downloading file {0}...', self.remote_file)
             if not self.local_file:
                 self.local_file = TempFile.get()
-            Utils.get_file_hash(self.remote_file, use_stored=True, save=True)
+            #Utils.get_file_hash(self.remote_file, use_stored=True, save=True)
             try:
                 Utils.copy_file(self.remote_file, self.local_file, retrials=1)
             except Utils.CopyException:
-                Utils.get_file_hash(self.remote_file, use_stored=False, save=True)
+             #   Utils.get_file_hash(self.remote_file, use_stored=False, save=True)
                 Utils.copy_file(self.remote_file, self.local_file, retrials=2)
 
             if self.data_convention == 'meteofrance':
diff --git a/earthdiagnostics/ocean/heatcontentlayer.py b/earthdiagnostics/ocean/heatcontentlayer.py
index 023909a09bee34676ad13ace072bedbbe1ae2bbf..bf0cc580a18e18a6377d81047ccd57d9d1e2b93e 100644
--- a/earthdiagnostics/ocean/heatcontentlayer.py
+++ b/earthdiagnostics/ocean/heatcontentlayer.py
@@ -4,10 +4,19 @@ import numpy as np
 
 from earthdiagnostics.box import Box
 from earthdiagnostics.constants import Basins
-from earthdiagnostics.diagnostic import Diagnostic, DiagnosticIntOption, DiagnosticBasinOption
+from earthdiagnostics.diagnostic import Diagnostic, DiagnosticIntOption, DiagnosticBasinListOption
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
+from bscearth.utils.log import Log
 
+import diagonals
+from diagonals import ohc
+from diagonals.mesh_helpers.nemo import Nemo
+
+import iris
+import iris.coords
+
+import netCDF4
 
 class HeatContentLayer(Diagnostic):
     """
@@ -36,15 +45,17 @@ class HeatContentLayer(Diagnostic):
     alias = 'ohclayer'
     "Diagnostic alias for the configuration file"
 
-    def __init__(self, data_manager, startdate, member, chunk, box, weight, min_level, max_level, data_convention):
+    def __init__(self, data_manager, startdate, member, chunk, box, areas,
+                 weight, layers, basins, data_convention):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
         self.member = member
         self.chunk = chunk
         self.box = box
+        self.areas = areas
         self.weight = weight
-        self.min_level = min_level
-        self.max_level = max_level
+        self.layers = layers
+        self.basins = basins
         self.required_vars = ['so', 'mlotst']
         self.generated_vars = ['scvertsum']
         self.data_convention = data_convention
@@ -65,62 +76,47 @@ class HeatContentLayer(Diagnostic):
         """
         options_available = (DiagnosticIntOption('min_depth'),
                              DiagnosticIntOption('max_depth'),
-                             DiagnosticBasinOption('basin', Basins().Global))
+                             DiagnosticBasinListOption('basins',
+                             Basins().Global))
         options = cls.process_options(options, options_available)
 
         box = Box(True)
         box.min_depth = options['min_depth']
         box.max_depth = options['max_depth']
+        basins = options['basins']
+        if not basins:
+            Log.error('Basins not recognized')
+            return ()
+
+        mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
+
+        layers = ((box.min_depth, box.max_depth),)
+
+        masks = {}
+        basins.sort()
+        for basin in basins:
+            masks[basin] = Utils.get_mask(basin)
+
+        areacello = mesh.get_areacello()
+        areas = ohc.get_basin_area(areacello, masks)
+
         job_list = list()
 
-        max_level, min_level, weight = cls._compute_weights(box)
+        e3t = mesh.get_k_length()
+        mask = mesh.get_landsea_mask()
+        depth = mesh.get_depth(cell_point='W')
+        weight = ohc.get_weights(layers, mask, e3t, depth)
+
+        del mask, depth, e3t
 
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job_list.append(HeatContentLayer(
                 diags.data_manager, startdate, member, chunk, box,
-                weight, min_level, max_level,
+                areas, weight, layers, basins,
                 diags.config.data_convention
             ))
         return job_list
 
-    @classmethod
-    def _compute_weights(cls, box):
-        depth, e3t, mask = cls._get_mesh_info()
-
-        def calculate_weight(e3t_point, depth_point, mask_point):
-            """Calculate the weight for each cell"""
-            if not mask_point:
-                return 0
-            top = depth_point
-            bottom = top + e3t_point
-            if bottom < box.min_depth or top > box.max_depth:
-                return 0
-            else:
-                if top < box.min_depth:
-                    top = box.min_depth
-                if bottom > box.max_depth:
-                    bottom = box.max_depth
-
-                return (bottom - top) * 1020 * 4000
-
-        calc = np.vectorize(calculate_weight, otypes='f')
-        weight = calc(e3t, depth, mask)
-        max_level, min_level = cls._get_used_levels(weight)
-        weight = weight[:, min_level:max_level, :]
-        return max_level, min_level, weight
-
-    @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.shape[1]
-        min_level = 0
-        while min_level < levels and not weight[:, min_level, :].any():
-            min_level += 1
-        max_level = min_level
-        while max_level < (levels - 1) and weight[:, max_level + 1, :].any():
-            max_level += 1
-        return max_level, min_level
-
     @classmethod
     def _get_mesh_info(cls):
         handler = Utils.open_cdf('mesh_zgr.nc')
@@ -150,42 +146,62 @@ class HeatContentLayer(Diagnostic):
 
     def request_data(self):
         """Request data required by the diagnostic"""
-        self.thetao = self.request_chunk(ModelingRealms.ocean, 'thetao', self.startdate, self.member, self.chunk)
+        self.thetao = self.request_chunk(ModelingRealms.ocean, 'thetao',
+                                         self.startdate, self.member,
+                                         self.chunk)
 
     def declare_data_generated(self):
         """Declare data to be generated by the diagnostic"""
-        self.heatc = self.declare_chunk(ModelingRealms.ocean, 'heatc', self.startdate, self.member, self.chunk,
-                                        box=self.box)
+        self.heatc = self.declare_chunk(ModelingRealms.ocean, 'heatc',
+                                        self.startdate, self.member,
+                                        self.chunk, box=self.box)
+        self.heatcsum = self.declare_chunk(ModelingRealms.ocean, 'heatcsum',
+                                           self.startdate, self.member,
+                                           self.chunk, box=self.box)
 
     def compute(self):
         """Run the diagnostic"""
-        nco = Utils.nco()
+        
         thetao_file = TempFile.get()
         results = TempFile.get()
-
+        results1D = TempFile.get()
         Utils.copy_file(self.thetao.local_file, thetao_file)
 
         handler = Utils.open_cdf(thetao_file)
         Utils.convert_units(handler.variables['thetao'], 'K')
-        heatc_sl = np.sum(handler.variables['thetao'][:, self.min_level:self.max_level, :] * self.weight, 1)
+        heatc_sl, heatc_sl1D = ohc.compute(self.layers, self.weight,
+                                           handler.variables['thetao'][:],
+                                           self.areas)
         handler.sync()
         handler.renameVariable('thetao', 'heatc_sl')
-        handler.close()
 
-        nco.ncks(
-            input=thetao_file,
-            output=results,
-            options=(
-                '-O -v {0.lon_name},{0.lat_name},time'.format(self.data_convention),
-            )
-        )
-        Utils.rename_variables(results, {'x': 'i', 'y': 'j'}, False)
-        handler_results = Utils.open_cdf(results)
-        var = handler_results.createVariable('heatc', float, ('time', 'j', 'i'), fill_value=1.e20)
+        results = TempFile.get()
+        handler_results = Utils.open_cdf(results, 'w')
+        Utils.copy_variable(handler, handler_results, 'time', True, True)
+        Utils.copy_variable(handler, handler_results, 'i', True, True)
+        Utils.copy_variable(handler, handler_results, 'j', True, True)
+#        Utils.rename_variables(results, {'x': 'i', 'y': 'j'}, False)
+        var = handler_results.createVariable('heatc', float,
+                                             ('time', 'j', 'i'),
+                                             fill_value=1.e20)
         var.units = 'J m-2'
         handler_results.sync()
-        handler_results.variables['heatc'][:] = heatc_sl
+        handler_results.variables['heatc'][:] = heatc_sl[0] # temporary fix, needs to loop over layers
         handler_results.close()
 
+        results1D = TempFile.get()
+        handler_results1D = Utils.open_cdf(results1D, 'w')
+        Utils.copy_variable(handler, handler_results1D, 'time', True, True)
+        handler_results1D.createDimension('region', len(self.basins))
+        handler_results1D.createDimension('region_length', 50)
+        var_region = handler_results1D.createVariable('region', 'S1', ('region', 'region_length'))
+        var_ohc1D = handler_results1D.createVariable('heatcsum', float, ('time', 'region',),)
+        handler_results1D.sync()
+        for i, basin in enumerate(self.basins):
+            var_region[i, ...] = netCDF4.stringtoarr(basin.name, 50)
+            var_ohc1D[..., i] = heatc_sl1D[i]
+        handler_results1D.close()
+
         Utils.setminmax(results, 'heatc')
         self.heatc.set_local_file(results)
+        self.heatcsum.set_local_file(results1D)
diff --git a/earthdiagnostics/ocean/moc.py b/earthdiagnostics/ocean/moc.py
index 7732726d661aca8aa4ed059dee409c54cd8af957..a294b824ce3a581533613de18d6e3b6d883dc9f6 100644
--- a/earthdiagnostics/ocean/moc.py
+++ b/earthdiagnostics/ocean/moc.py
@@ -67,8 +67,10 @@ class Moc(Diagnostic):
     def __eq__(self, other):
         if self._different_type(other):
             return False
-        return self.startdate == other.startdate and \
-               self.member == other.member and self.chunk == other.chunk
+        return (
+            self.startdate == other.startdate and
+            self.member == other.member and self.chunk == other.chunk
+        )
 
     @classmethod
     def generate_jobs(cls, diags, options):
@@ -136,7 +138,6 @@ class Moc(Diagnostic):
         del vo, e1v, e3v
         self._save_result(moc_results, mesh)
 
-
     def _save_result(self, result, mesh):
         temp = TempFile.get()
         handler_source = Utils.open_cdf(self.variable_file.local_file)
@@ -151,7 +152,6 @@ class Moc(Diagnostic):
         handler_temp.createDimension('region', len(result))
         handler_temp.createDimension('region_length', 50)
 
-
         var_region = handler_temp.createVariable('region', 'S1',
                                                  ('region', 'region_length'))
 
diff --git a/earthdiagnostics/ocean/regionmean.py b/earthdiagnostics/ocean/regionmean.py
index 792793a535e7f2b9c0d0743ee206541900717677..89a14fbcf5da07a944e90b71c05b04da88a75f9f 100644
--- a/earthdiagnostics/ocean/regionmean.py
+++ b/earthdiagnostics/ocean/regionmean.py
@@ -22,7 +22,6 @@ import diagonals.regmean as regmean
 from diagonals.mesh_helpers.nemo import Nemo
 
 
-
 class RegionMean(Diagnostic):
     """
     Computes the mean value of the field (3D, weighted).
@@ -79,10 +78,10 @@ class RegionMean(Diagnostic):
             self.box == other.box and self.variable == other.variable
 
     def __str__(self):
-        return 'Region mean Startdate: {0.startdate} Member: {0.member}' \
-               'Chunk: {0.chunk} Variable: {0.variable} ' \
-               'Box: {0.box} Save 3D: {0.save3d} Save variance: {0.variance}' \
-               'Grid point: {0.grid_point}'.format(self)
+        return ('Region mean Startdate: {0.startdate} Member: {0.member} '
+                'Chunk: {0.chunk} Variable: {0.variable} '
+                'Box: {0.box} Save 3D: {0.save3d} Save variance: {0.variance} '
+                'Grid point: {0.grid_point}'.format(self))
 
     def __hash__(self):
         return hash(str(self))
@@ -101,8 +100,7 @@ class RegionMean(Diagnostic):
         options_available = (DiagnosticDomainOption(),
                              DiagnosticVariableOption(diags.data_manager.config.var_manager),
                              DiagnosticOption('grid_point', 'T'),
-                             DiagnosticBasinListOption('basins',
-                                                       Basins().Global),
+                             DiagnosticBasinListOption('basins', 'global'),
                              DiagnosticIntOption('min_depth', -1),
                              DiagnosticIntOption('max_depth', -1),
                              DiagnosticIntOption('min_lat', -1),
@@ -117,17 +115,18 @@ class RegionMean(Diagnostic):
         box = Box()
         box.min_depth = options['min_depth']
         box.max_depth = options['max_depth']
-        box.min_lat = options['min_lat']
-        box.max_lat = options['max_lat']
-        box.min_lon = options['min_lon']
-        box.max_lon = options['max_lon']
+        if options['min_lat'] != -1:
+            box.min_lat = options['min_lat']
+            box.max_lat = options['max_lat']
+        if options['min_lon'] != -1 or options['max_lon'] != -1:
+            box.min_lon = options['min_lon']
+            box.max_lon = options['max_lon']
 
         basins = options['basins']
         if not basins:
             Log.error('Basins not recognized')
             return()
 
-
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job = RegionMean(diags.data_manager, startdate, member, chunk,
@@ -165,10 +164,13 @@ class RegionMean(Diagnostic):
             masks[basin] = Utils.get_mask(basin)
 
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
-        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and \
-            self.box.min_lon is not -1 and self.box.max_lat is not -1:
-            name = '{0}_{1}'.format(Box.get_lat_str(self.box),
-                                    Box.get_lon_str(self.box))
+        if (
+            self.box.min_lat is not -1 and self.box.max_lat is not -1 and
+            self.box.min_lon is not -1 and self.box.max_lat is not -1
+        ):
+            name = '{0}_{1}'.format(
+                Box.get_lat_str(self.box), Box.get_lon_str(self.box)
+            )
 
             masks[name] = mesh.get_region_mask(self.box.min_lat,
                                                self.box.max_lat,
@@ -184,7 +186,6 @@ class RegionMean(Diagnostic):
         mean = regmean.compute_regmean_2D(data.data, masks, areacello)
         self._save_result_2D('mean', mean, data)
 
-
     def _meand_3d_variable(self, data, mesh, masks):
         areacello = mesh.get_areacello(cell_point=self.grid_point)
         e3 = self._try_load_cube(3)
@@ -279,7 +280,6 @@ class RegionMean(Diagnostic):
                                                        box=box_save,
                                                        region=self.basins)
 
-
     def _save_result_2D(self, var, result, data):
         final_name = '{1}{0}'.format(var, self.variable)
         temp = TempFile.get()
@@ -299,7 +299,6 @@ class RegionMean(Diagnostic):
         handler_temp.close()
         self.declared[final_name].set_local_file(temp, diagnostic=self)
 
-
     def _save_result_3D(self, var, result, data):
         final_name = '{1}3d{0}'.format(var, self.variable)
         temp = TempFile.get()
@@ -321,11 +320,10 @@ class RegionMean(Diagnostic):
 
         var = handler_temp.createVariable('{1}3d{0}'.format(var, self.variable),
                                           float, ('time', 'lev', 'region',),)
-                                          
+
         var.units = '{0}'.format(data.units)
         for i, basin in enumerate(result):
             var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
             var[..., i] = result[basin]
         handler_temp.close()
         self.declared[final_name].set_local_file(temp, diagnostic=self)
-
diff --git a/earthdiagnostics/ocean/regionsum.py b/earthdiagnostics/ocean/regionsum.py
index f9207e5999dcf6c015ac15c927ab43d366de8bb9..99b46ce19a28bccf74372b509098359e62a09776 100644
--- a/earthdiagnostics/ocean/regionsum.py
+++ b/earthdiagnostics/ocean/regionsum.py
@@ -24,6 +24,7 @@ from earthdiagnostics.utils import Utils, TempFile
 import diagonals.regsum as regsum
 from diagonals.mesh_helpers.nemo import Nemo
 
+
 class RegionSum(Diagnostic):
     """
     Computes the sum of the field (3D, weighted).
@@ -88,11 +89,9 @@ class RegionSum(Diagnostic):
                'Grid point: {0.grid_point} Box: {0.box} Save 3D: {0.save3d}' \
                'Original grid: {0.grid} Basin: {0.basins}'.format(self)
 
-
     def __hash__(self):
         return hash(str(self))
 
-
     @classmethod
     def generate_jobs(cls, diags, options):
         """
@@ -167,7 +166,8 @@ class RegionSum(Diagnostic):
             masks[basin] = Utils.get_mask(basin)
 
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
-        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and self.box.min_lon is not -1 and self.box.max_lat is not -1:
+        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and \
+           self.box.min_lon is not -1 and self.box.max_lat is not -1:
             name = '{0}_{1}'.format(Box.get_lat_str(self.box), Box.get_lon_str(self.box))
 
             masks[name] = mesh.get_region_mask(self.box.min_lat,
@@ -179,13 +179,11 @@ class RegionSum(Diagnostic):
         else:
             self._sum_2d_var(data, mesh, masks)
 
-
     def _sum_2d_var(self, data, mesh, masks):
         areacello = mesh.get_areacello(cell_point=self.grid_point)
         varsum = regsum.compute_regsum_2D(data.data, masks, areacello)
         self._save_result_2D('sum', varsum, data)
 
-
     def _sum_3d_variable(self, data, mesh, masks):
         areacello = mesh.get_areacello(cell_point=self.grid_point)
         tmask = iris.load_cube('mesh_hgr.nc', '{0}mask'.format(self.grid_point))
@@ -206,7 +204,6 @@ class RegionSum(Diagnostic):
                                                     volcello, tmask)
             self._save_result_3D('sum', varsum3d, data)
 
-
     def _try_load_cube(self, number):
         try:
             cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}'.format(number, self.grid_point))
@@ -239,7 +236,6 @@ class RegionSum(Diagnostic):
         data = iris.load_cube(self.variable_file.local_file)
         return self._rename_depth(data)
 
-
     def _rename_depth(self, data):
         for coord_name in ('model_level_number', 'Vertical T levels', 'lev'):
             if data.coords(coord_name):
@@ -249,7 +245,6 @@ class RegionSum(Diagnostic):
                 break
         return data
 
-
     def _fix_file_metadata(self):
         handler = Utils.open_cdf(self.variable_file.local_file)
         var = handler.variables[self.variable]
@@ -314,7 +309,6 @@ class RegionSum(Diagnostic):
         self.declared[final_name] = self.declare_chunk(ModelingRealms.ocean, final_name, self.startdate, self.member,
                                                        self.chunk, box=box_save, region=self.basins, grid=self.grid)
 
-
     def _save_result_2D(self, var, result, data):
         final_name = '{1}{0}'.format(var, self.variable)
         temp = TempFile.get()
@@ -334,7 +328,6 @@ class RegionSum(Diagnostic):
         handler_temp.close()
         self.declared[final_name].set_local_file(temp, diagnostic=self)
 
-
     def _save_result_3D(self, var, result, data):
         final_name = '{1}3d{0}'.format(var, self.variable)
         temp = TempFile.get()
diff --git a/earthdiagnostics/ocean/siasiesiv.py b/earthdiagnostics/ocean/siasiesiv.py
index 5cf4a9050cf6e4fa6c905a20d953d7e387d4ccbe..94c65d3c98dcc1039169474af7f6f0c8dc034904 100644
--- a/earthdiagnostics/ocean/siasiesiv.py
+++ b/earthdiagnostics/ocean/siasiesiv.py
@@ -70,6 +70,9 @@ class Siasiesiv(Diagnostic):
                'Basins: {1} Omit volume: {0.omit_volume}'.format(self,
                                                                  ','.join(str(basin) for basin in self.masks.keys()))
 
+    def __hash__(self):
+        return hash(str(self))
+
     @classmethod
     def generate_jobs(cls, diags, options):
         """
@@ -140,8 +143,10 @@ class Siasiesiv(Diagnostic):
             self.sic.local_file, self.sic_varname
         )
         sic = iris.load_cube(self.sic.local_file)
-        if sic.units.origin == '%' and sic.data.max() < 2:
-            sic.units = '1.0'
+        #if sic.units.origin == '%' and sic.data.max() < 2:
+        #    sic.units = '1.0'
+        if sic.units.origin is not '%': 
+            sic.convert_units('%')
 
         sic_slices = []
         for sic_data in sic.slices_over('time'):
@@ -150,7 +155,17 @@ class Siasiesiv(Diagnostic):
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
         areacello = mesh.get_areacello(cell_point='T')
         gphit = mesh.get_grid_latitude(cell_point='T')
-        self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'] = siasie.compute(gphit, areacello, sic_slices, self.masks)
+
+        if not self.omit_volume:
+            sit = iris.load_cube(self.sit.local_file)
+            sit_slices = []
+            for sit_data in sit.slices_over('time'):
+                sit_data.data = np.ma.filled(sit_data.data, 0.0).astype(np.float32)
+                sit_slices.append(sit_data)
+            self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'], self.results['sivoln'], self.results['sivols'] = siasie.compute(gphit, areacello, sic_slices, self.masks, sit_slices)
+        else:
+            self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'] = siasie.compute(gphit, areacello, sic_slices, self.masks, None)
+
 
         self.save()
 
@@ -179,7 +194,10 @@ class Siasiesiv(Diagnostic):
                                                  ('region', 'region_length'))
             var_res = handler_temp.createVariable('{0}'.format(var), float,
                                                   ('time', 'region',))
-            var_res.units = 'm^2'
+            if var in ('sivoln', 'sivols'):
+                var_res.units = 'm^3'
+            else:
+                var_res.units = 'm^2'
             for i, basin in enumerate(self.masks):
                 var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
                 var_res[..., i] = res[i, ...]
diff --git a/earthdiagnostics/ocean/zonalmean.py b/earthdiagnostics/ocean/zonalmean.py
index 5e9177907d36fcf050f8a36e8693e8a5ac335a27..b7b9908e970ed95c25c71071a23b36dd292df3ba 100644
--- a/earthdiagnostics/ocean/zonalmean.py
+++ b/earthdiagnostics/ocean/zonalmean.py
@@ -11,13 +11,18 @@ from iris.cube import Cube, CubeList
 import numpy as np
 import numba
 
+import netCDF4
+
 from earthdiagnostics.box import Box
 from earthdiagnostics.constants import Basins
 from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, DiagnosticIntOption, DiagnosticDomainOption, \
-    DiagnosticBoolOption, DiagnosticBasinOption, DiagnosticVariableOption
+    DiagnosticBoolOption, DiagnosticBasinListOption, DiagnosticVariableOption
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
 
+import diagonals.zonmean as zonmean
+from diagonals.mesh_helpers.nemo import Nemo
+
 
 class ZonalMean(Diagnostic):
     """
@@ -40,14 +45,14 @@ class ZonalMean(Diagnostic):
     alias = 'zonmean'
     "Diagnostic alias for the configuration file"
 
-    def __init__(self, data_manager, startdate, member, chunk, domain, variable, basin, grid_point):
+    def __init__(self, data_manager, startdate, member, chunk, domain, variable, basins, grid_point):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
         self.member = member
         self.chunk = chunk
         self.domain = domain
         self.variable = variable
-        self.basin = basin
+        self.basins = basins
         self.grid_point = grid_point
 
         self.declared = {}
@@ -83,14 +88,14 @@ class ZonalMean(Diagnostic):
             DiagnosticDomainOption(),
             DiagnosticVariableOption(diags.data_manager.config.var_manager),
             DiagnosticOption('grid_point', 'T'),
-            DiagnosticBasinOption('basin', Basins().Global),
+            DiagnosticBasinListOption('basins', Basins().Global),
         )
         options = cls.process_options(options, options_available)
 
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job = ZonalMean(diags.data_manager, startdate, member, chunk,
-                            options['domain'], options['variable'], options['basin'],
+                            options['domain'], options['variable'], options['basins'],
                             options['grid_point'].lower())
             job_list.append(job)
 
@@ -106,76 +111,24 @@ class ZonalMean(Diagnostic):
         """Declare data to be generated by the diagnostic"""
         self.zonal_mean = self.declare_chunk(
             ModelingRealms.ocean, self.variable + 'zonal',
-            self.startdate, self.member, self.chunk, region=self.basin
+            self.startdate, self.member, self.chunk, region=self.basins
         )
 
     def compute(self):
         """Run the diagnostic"""
         self._fix_file_metadata()
+        mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
+        areacello = mesh.get_areacello(cell_point=self.grid_point)
+        lats = mesh.get_mesh_var('latitude', dtype=np.float32)
         data = self._load_data()
-        self._meand_3d_variable(data)
-
-    def _meand_3d_variable(self, data):
-        e1 = self._try_load_cube(1)
-        e2 = self._try_load_cube(2)
-        mask = np.squeeze(Utils.get_mask(self.basin, True))
-        mask = e1.data * e2.data * mask
-        if len(mask.shape) == 2:
-            data.add_aux_coord(
-                iris.coords.AuxCoord(mask.data, long_name='mask'),
-                data.coord_dims('latitude')
-            )
-        else:
-            data.add_aux_coord(
-                iris.coords.AuxCoord(mask.data, long_name='mask'),
-                data.coord_dims('depth') + data.coord_dims('latitude')
-            )
-
-        @numba.njit()
-        def get_zonal_mean(variable, weight, latitude):
-            total = np.zeros(180, np.float64)
-            weights = np.zeros(180, np.float64)
-            for i in range(variable.shape[0]):
-                for j in range(variable.shape[1]):
-                    if weight[i, j] == 0:
-                        continue
-                    bin_value = int(round(latitude[i, j]) + 90)
-                    weights[bin_value] += weight[i, j]
-                    total[bin_value] += variable[i, j] * weight[i, j]
-            return total / weights
-
-        mean = iris.cube.CubeList()
-        lat_coord = None
-        for map_slice in data.slices_over('time'):
-            # Force data loading
-            map_slice.data
-            surface_cubes = iris.cube.CubeList()
-            for surface_slice in map_slice.slices_over('depth'):
-                value = get_zonal_mean(
-                    surface_slice.data,
-                    surface_slice.coord('mask').points,
-                    surface_slice.coord('latitude').points,
-                )
-                cube = Cube(value)
-                cube.add_aux_coord(surface_slice.coord('depth'))
-                if lat_coord is None:
-                    lat_coord = surface_slice.coord('latitude')
-                    lat_coord = lat_coord.copy(
-                        np.arange(-90, 90, dtype=np.float32)
-                    )
-                    lat_coord = iris.coords.DimCoord.from_coord(lat_coord)
-                cube.add_dim_coord(lat_coord, 0)
-                surface_cubes.append(cube)
-            time_cube = surface_cubes.merge_cube()
-            time_cube.add_aux_coord(map_slice.coord('time'))
-            mean.append(time_cube)
-        cube = mean.merge_cube()
-        cube.var_name = 'result'
-        cube.units = data.units
-        cube.attributes = data.attributes
-        temp = TempFile.get()
-        iris.save(cube, temp)
-        self.zonal_mean.set_local_file(temp, rename_var='result', region=self.basin)
+        masks = {}
+        self.basins.sort()
+        for basin in self.basins:
+            masks[basin] = Utils.get_mask(basin)
+        area_basin = zonmean.get_basin_area(areacello, masks)
+        zonalmean = zonmean.compute_zonmean(data.data, lats, area_basin)
+        self._save_result(zonalmean, data)
+
 
     def _try_load_cube(self, number):
         try:
@@ -231,3 +184,40 @@ class ZonalMean(Diagnostic):
         var.coordinates = coordinates
         handler.close()
         return has_levels
+
+
+    def _save_result(self, result, data):
+        final_name = '{0}zonal'.format(self.variable)
+        temp = TempFile.get()
+        handler_source = Utils.open_cdf(self.variable_file.local_file)
+        handler_temp = Utils.open_cdf(temp, 'w')
+        Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
+        handler_temp.createDimension('region', len(result))
+        handler_temp.createDimension('region_length', 50)
+        handler_temp.createDimension('lev', data.shape[1])
+        handler_temp.createDimension('lat', 180)
+
+        var_lat = handler_temp.createVariable('lat', float, 'lat')
+        var_lat[...] = np.arange(-90, 90, dtype=np.float32)
+        var_lat.units = 'degrees_north'
+        var_lat.long_name = 'latitude'
+        var_lat.standard_name = 'latitude'
+
+        var_level = handler_temp.createVariable('lev', float, 'lev')
+        var_level[...] = data.coord('depth').points
+        var_level.units = 'm'
+        var_level.axis = 'Z'
+        var_level.positive = 'down'
+        var_level.long_name = 'ocean depth coordinate'
+        var_level.standard_name = 'depth'
+        var_region = handler_temp.createVariable('region', 'S1',
+                                                 ('region', 'region_length'))
+        var = handler_temp.createVariable('{0}zonal'.format(self.variable),
+                                          float, ('time', 'lev', 'lat','region',),)
+        var.units = '{0}'.format(data.units)
+        for i, basin in enumerate(result):
+            var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
+            var[..., i] = result[basin]
+        handler_temp.close()
+        self.declared[final_name].set_local_file(temp, diagnostic=self)
+
diff --git a/earthdiagnostics/utils.py b/earthdiagnostics/utils.py
index c08c1f08b84e1cf27fbcdce11e6e8dc3bde99274..0920a3657014f9d13d2af41754ff39e6537d524d 100644
--- a/earthdiagnostics/utils.py
+++ b/earthdiagnostics/utils.py
@@ -328,21 +328,21 @@ class Utils(object):
                 # This can be due to a race condition. If directory already exists, we don have to do nothing
                 if not os.path.exists(dirname_path):
                     raise ex
-        hash_destiny = None
-        Log.debug('Hashing original file... {0}', datetime.datetime.now())
-        hash_original = Utils.get_file_hash(source, use_stored=use_stored_hash)
-
-        if retrials < 1:
-            retrials = 1
-
-        while hash_original != hash_destiny:
-            if retrials == 0:
-                raise Utils.CopyException('Can not copy {0} to {1}'.format(source, destiny))
-            Log.debug('Copying... {0}', datetime.datetime.now())
-            shutil.copyfile(source, destiny)
-            Log.debug('Hashing copy ... {0}', datetime.datetime.now())
-            hash_destiny = Utils.get_file_hash(destiny, save=save_hash)
-            retrials -= 1
+        #hash_destiny = None
+        #Log.debug('Hashing original file... {0}', datetime.datetime.now())
+        #hash_original = Utils.get_file_hash(source, use_stored=use_stored_hash)
+
+       # if retrials < 1:
+       #     retrials = 1
+
+       # while hash_original != hash_destiny:
+       #     if retrials == 0:
+       #         raise Utils.CopyException('Can not copy {0} to {1}'.format(source, destiny))
+        Log.debug('Copying... {0}', datetime.datetime.now())
+        shutil.copyfile(source, destiny)
+       #     Log.debug('Hashing copy ... {0}', datetime.datetime.now())
+            #hash_destiny = Utils.get_file_hash(destiny, save=save_hash)
+        #    retrials -= 1
         Log.debug('Finished {0}', datetime.datetime.now())
 
     @staticmethod
@@ -426,8 +426,8 @@ class Utils(object):
         time.sleep(2)
         shutil.rmtree(source)
 
-    @staticmethod
-    def get_file_hash(filepath, use_stored=False, save=False):
+ #   @staticmethod
+   # def get_file_hash(filepath, use_stored=False, save=False):
         """
         Get the xxHash hash for a given file
 
@@ -439,33 +439,33 @@ class Utils(object):
         save: bool, optional
             If True, saves the hash to a file
         """
-        if use_stored:
-            hash_file = Utils._get_hash_filename(filepath)
-            if os.path.isfile(hash_file):
-                hash_value = open(hash_file, 'r').readline()
-                return hash_value
-
-        blocksize = 104857600
-        hasher = xxhash.xxh64()
-        with open(filepath, 'rb') as afile:
-            buf = afile.read(blocksize)
-            while len(buf) > 0:
-                hasher.update(buf)
-                buf = afile.read(blocksize)
-        hash_value = hasher.hexdigest()
-        if save:
-            hash_file = open(Utils._get_hash_filename(filepath), 'w')
-            hash_file.write(hash_value)
-            hash_file.close()
-
-        return hash_value
-
-    @staticmethod
-    def _get_hash_filename(filepath):
-        folder = os.path.dirname(filepath)
-        filename = os.path.basename(filepath)
-        hash_file = os.path.join(folder, '.{0}.xxhash64.hash'.format(filename))
-        return hash_file
+      #  if use_stored:
+      #      hash_file = Utils._get_hash_filename(filepath)
+      #      if os.path.isfile(hash_file):
+      #          hash_value = open(hash_file, 'r').readline()
+      #          return hash_value
+
+     #   blocksize = 104857600
+     #   hasher = xxhash.xxh64()
+     #   with open(filepath, 'rb') as afile:
+     #       buf = afile.read(blocksize)
+     #       while len(buf) > 0:
+     #           hasher.update(buf)
+     #           buf = afile.read(blocksize)
+     #   hash_value = hasher.hexdigest()
+     #   if save:
+     #       hash_file = open(Utils._get_hash_filename(filepath), 'w')
+     #       hash_file.write(hash_value)
+     #       hash_file.close()
+
+    #    return hash_value
+
+  #  @staticmethod
+  #  def _get_hash_filename(filepath):
+  #      folder = os.path.dirname(filepath)
+  #      filename = os.path.basename(filepath)
+  #      hash_file = os.path.join(folder, '.{0}.xxhash64.hash'.format(filename))
+  #      return hash_file
 
     @staticmethod
     def execute_shell_command(command, log_level=Log.DEBUG):
diff --git a/test/unit/data_convention/test_primavera.py b/test/unit/data_convention/test_primavera.py
index f2e6739b04288f629c54d64ca6fed90b12dd59b8..52a06faadfdec95a4e5d6b817631ee7fd76521a0 100644
--- a/test/unit/data_convention/test_primavera.py
+++ b/test/unit/data_convention/test_primavera.py
@@ -32,6 +32,7 @@ class TestPrimaveraConvention(TestCase):
 
         self.config.cmor.activity = 'activity'
         self.config.cmor.append_startdate = False
+        self.config.cmor.append_startdate_year_only = False
         self.config.cmor.initialization_number = 1
         self.config.cmor.version = 'version'
         self.config.cmor.default_ocean_grid = 'ocean_grid'
@@ -148,6 +149,19 @@ class TestPrimaveraConvention(TestCase):
         self.assertEqual(file_path,
                          'var_Omon_model_experiment_name_s19900101-r2i1p1f1_ocean_grid_199001-199001.nc')
 
+    def test_get_filename_with_startdate_year_only(self):
+        """Test get_filename with startdate"""
+        self.config.cmor.append_startdate = True
+        self.config.cmor.append_startdate_year_only = True
+        cmor_var = Mock()
+        omon = Mock()
+        omon.name = 'Omon'
+        cmor_var.get_table.return_value = omon
+        file_path = self.convention.get_file_name('19900101', 1, ModelingRealms.ocean, 'var', cmor_var,
+                                                  Frequencies.monthly, 1, None, None, None)
+        self.assertEqual(file_path,
+                         'var_Omon_model_experiment_name_s1990-r2i1p1f1_ocean_grid_199001-199001.nc')
+
     def test_get_filename_no_cmor_var(self):
         """Test get_filename not passing cmor_var"""
         cmor_var = Mock()
diff --git a/test/unit/ocean/test_region_mean.py b/test/unit/ocean/test_region_mean.py
index 966ed536c97807961da450d383dad07558c66a30..35f8a09cb38d5e64793fac60d43e1c912409b9ae 100644
--- a/test/unit/ocean/test_region_mean.py
+++ b/test/unit/ocean/test_region_mean.py
@@ -38,23 +38,23 @@ class TestRegionMean(TestCase):
         jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var'])
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 't'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 't'))
 
         jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U'])
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
 
         jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global'])
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
 
         box = Box()
         box.min_depth = 1.0
@@ -63,27 +63,34 @@ class TestRegionMean(TestCase):
         jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10'])
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
 
-        jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', 'false'])
+        jobs = RegionMean.generate_jobs(
+            self.diags,
+            ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', '', '', '', '', 'false']
+        )
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, False, Basins().Global, False, ''))
+                                             box, False, False, Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, False, Basins().Global, False, ''))
+                                             box, False, False, Basins().Global, 'u'))
 
-        jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', 'false',
-                                                     'True'])
+        jobs = RegionMean.generate_jobs(
+            self.diags,
+            ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', '', '', '', '', 'false', 'True']
+        )
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, False, Basins().Global, True, ''))
+                                             box, False, True,  Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, False, Basins().Global, True, ''))
+                                             box, False, True, Basins().Global, 'u'))
 
-        jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', 'false',
-                                                     'True', 'grid'])
+        jobs = RegionMean.generate_jobs(
+            self.diags,
+            ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', '', '', '', '', 'false', 'True', 'grid']
+        )
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
                                              box, False, Basins().Global, True, 'grid'))
@@ -94,8 +101,11 @@ class TestRegionMean(TestCase):
             RegionMean.generate_jobs(self.diags, ['diagnostic'])
 
         with self.assertRaises(DiagnosticOptionError):
-            RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', 'false',
-                                                  'True', 'grid', 'extra'])
+            RegionMean.generate_jobs(
+                self.diags,
+                ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', '', '', '', '', 'false',
+                 'True', 'grid', 'extra']
+            )
 
     def test_str(self):
         box = Box()