diff --git a/VERSION b/VERSION
index fa8da2074bb62a20ddd742cf2e020bb750e75e78..62841113a966235da11a00d74683a1df0754aed7 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-3.5.8
+3.5.9
diff --git a/doc/source/conf.py b/doc/source/conf.py
index cafff619ea96be1d926b16db1ac2aa70133d09f6..b450227a83c51fd0390f9962c1b243527c7e5fe9 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -65,7 +65,7 @@ copyright = u"2020, BSC-CNS Earth Sciences Department"
 # The short X.Y version.
 version = "3.5"
 # The full version, including alpha/beta/rc tags.
-release = "3.5.8"
+release = "3.5.9"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/earthdiagnostics/ocean/interpolatecdo.py b/earthdiagnostics/ocean/interpolatecdo.py
index 9873f8c8ba47b0df4a3c61acf49d633b692cfad9..12421854d2c22ae81912b0e5de63f813cf21aa8c 100644
--- a/earthdiagnostics/ocean/interpolatecdo.py
+++ b/earthdiagnostics/ocean/interpolatecdo.py
@@ -305,6 +305,8 @@ class InterpolateCDO(Diagnostic):
         except AttributeError:
             units = None
         coordinates = list()
+        if 'type' in var.coordinates:
+            coordinates.append('type')
         for dim in var.dimensions:
             if dim == "i":
                 coordinates.append(lon_name)
diff --git a/earthdiagnostics/ocean/moc.py b/earthdiagnostics/ocean/moc.py
index 4f3f3c5dd5789659b162214d1c1f26bda0037b67..1d47366363653a58b1c4ef6f3152895e94cd61a3 100644
--- a/earthdiagnostics/ocean/moc.py
+++ b/earthdiagnostics/ocean/moc.py
@@ -149,11 +149,10 @@ class Moc(Diagnostic):
         handler_source = Utils.open_cdf(self.variable_file.local_file)
         handler_temp = Utils.open_cdf(temp, "w")
         gphiv = np.squeeze(mesh.get_grid_latitude(cell_point="V"))
-        max_gphiv = np.unravel_index(np.argmax(gphiv), gphiv.shape)[1]
+        max_gphiv = np.unravel_index(np.argmax(gphiv), gphiv.shape)[1] - 1
 
         Utils.copy_variable(handler_source, handler_temp, "time", True, True)
         Utils.copy_variable(handler_source, handler_temp, "lev", True, True)
-        handler_temp.createDimension("i", 1)
         handler_temp.createDimension("j", gphiv.shape[0])
         handler_temp.createDimension("region", len(result))
         handler_temp.createDimension("region_length", 50)
@@ -163,24 +162,18 @@ class Moc(Diagnostic):
         )
 
         lat = handler_temp.createVariable(
-            self.data_convention.lat_name, float, ("j", "i"))
+            self.data_convention.lat_name, float, ("j", ))
         lat[...] = gphiv[:, max_gphiv]
         lat.units = "degrees_north"
         lat.long_name = "Latitude"
         lat.standard_name = "latitude"
-
-        lon = handler_temp.createVariable(
-            self.data_convention.lon_name, float, ("j", "i"))
-        lon[...] = 0
-        lon.units = "degrees_east"
-        lon.long_name = "Longitude"
-        lon.standard_name = "longitude"
+        lat.axis = 'Y'
 
         var = handler_temp.createVariable(
-            "vsftmyz", float, ("time", "lev", "i", "j", "region")
+            "vsftmyz", float, ("time", "lev", "j", "region")
         )
         var.units = "Sverdrup"
-        var.coordinates = "lev time latitude longitude"
+        var.coordinates = "latitude region"
         var.long_name = "Ocean meridional overturning volume streamfunction"
         var.missing_value = 1e20
         var.fill_value = 1e20