diff --git a/Makefile b/Makefile
index 7c493d8d12ad97b7b73919a3bd8d01c97735e052..0578d4bafa75b11b101807d64b324628bdc82309 100644
--- a/Makefile
+++ b/Makefile
@@ -94,13 +94,14 @@ full_installation_for_dev: create_env set_env_vars_for_dev install_software_for_
 	Run 'conda activate $(CONDA_ENV_PATH)' to use the environment."
 
 # Combined task to make full local installation for user
-full_local_installation: create_env install_software run_tests
+full_local_installation: create_conda_env install_software run_tests
 	@echo "Local installation and testing completed. \
 	Run 'conda activate $(CONDA_ENV_PATH)' to use the environment."
 
 # Task to clean up the environment
 clean:
 	conda env remove -p $(CONDA_ENV_PATH) -y
+	conda clean --all --yes
 
 all: full_installation
 	@echo "Installation completed. Run 'conda activate $(CONDA_ENV_PATH)' to use the environment.\
diff --git a/environment.yml b/environment.yml
index 788450fd74b9612639bf30e7c046c49e35d7eb69..c49e204f5e011c000b00e02e817abe547bbd72d0 100755
--- a/environment.yml
+++ b/environment.yml
@@ -7,7 +7,6 @@ dependencies:
   - libnetcdf=*=mpi_mpich*
   - netCDF4=*=mpi_mpich*
   - h5py=*=mpi_mpich*
-  - pytest
   - pytest-cov
   - pycodestyle>=2.10.0
   - geopandas>=0.10.2
@@ -16,7 +15,8 @@ dependencies:
   - pyproj~=3.2.1
   - setuptools>=66.1.1
   - pytest>=7.2.1
-  - shapely
+  - shapely<2.0
+  - pygeos>=0.10.0
   - mpi4py ~= 3.1.4
   - eccodes
   - python-eccodes
diff --git a/nes/methods/horizontal_interpolation.py b/nes/methods/horizontal_interpolation.py
index 36a16b510cdfd0ed8e26785d47de820c09bb598e..717a3b0c2e013bc039a5e833e92daf040ec5ec2b 100644
--- a/nes/methods/horizontal_interpolation.py
+++ b/nes/methods/horizontal_interpolation.py
@@ -58,7 +58,7 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind="Neares
         weights, idx = __get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours,
                                                 only_create_wm, wm, flux)
     elif self.parallel_method in ["Y", "X"]:
-        print(weight_matrix_path)
+
         weights, idx = __get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours,
                                                  only_create_wm, wm, flux)
     else:
@@ -347,7 +347,6 @@ def __get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighb
         weight_matrix = wm
 
     elif weight_matrix_path is not None:
-        print(weight_matrix)
         with FileLock(weight_matrix_path + "{0:03d}.lock".format(self.rank)):
             if os.path.isfile(weight_matrix_path):
                 if self.master:
diff --git a/nes/methods/spatial_join.py b/nes/methods/spatial_join.py
index cc96975e2c6b7dd6e1d1feade6bf398d74e3d3f1..bfe9cc16ad09b09091c1958c0d99ac9fc25569b8 100644
--- a/nes/methods/spatial_join.py
+++ b/nes/methods/spatial_join.py
@@ -89,9 +89,9 @@ def __prepare_external_shapefile(self, ext_shp, var_list, info=False, apply_bbox
     GeoDataFrame
         External shapefile.
     """
-    if var_list is None:
-        var_list = []
-    
+    if isinstance(var_list, str):
+        var_list = [var_list]
+
     if isinstance(ext_shp, str):
         # Reading external shapefile
         if self.master and info:
@@ -101,7 +101,8 @@ def __prepare_external_shapefile(self, ext_shp, var_list, info=False, apply_bbox
         else:
             ext_shp = read_file(ext_shp)
 
-        ext_shp = ext_shp[var_list + ["geometry"]]
+        if var_list is not None:
+            ext_shp = ext_shp[var_list + ["geometry"]]
     else:
         msg = "WARNING!!! "
         msg += "External shapefile already read. If you pass the path to the shapefile instead of the opened shapefile "
diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py
index f97eacaea7945e3d613bd6a8f6f37319339180fb..59b27fb81f911d22f134792d7c9af7198bec3a46 100644
--- a/nes/nc_projections/default_nes.py
+++ b/nes/nc_projections/default_nes.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 
-import os, sys
+import os
+import sys
 from gc import collect
 from warnings import warn
 from math import isclose
@@ -1407,39 +1408,155 @@ class Nes(object):
 
         return None
 
+    def get_global_limits(self, limits: dict) -> dict:
+        """
+        Compute the global axis limits across all MPI ranks.
+
+        This function uses MPI reduce operations to determine the minimum and maximum
+        values for each axis across all ranks. `None` values are temporarily replaced
+        with `inf` or `-inf` to avoid issues during reduction and are restored afterward.
+
+        Parameters
+        ----------
+        limits : dict
+            A dictionary containing local axis limits for the current MPI rank.
+            Expected keys:
+            - `x_min`, `x_max`
+            - `y_min`, `y_max`
+            - `z_min`, `z_max`
+            - `t_min`, `t_max`
+            Values can be numerical or `None` (which will be processed accordingly).
+
+        Returns
+        -------
+        dict
+            A dictionary with the global minimum and maximum values for each axis,
+            computed across all MPI ranks. `None` is returned for dimensions that
+            were originally `None`.
+
+        Notes
+        -----
+        - The function assumes `self.comm` is a valid `MPI.COMM_WORLD` communicator.
+        - The reduction is performed only for numerical values; `None` values are handled explicitly.
+        - The final dictionary is only meaningful in rank 0. Other ranks return an undefined result.
+
+        Examples
+        --------
+        Suppose we have four MPI ranks with the following `limits`:
+
+        Rank 0:
+            {'x_min': -91, 'x_max': 50, 'y_min': 450, 'y_max': 475, 'z_min': None, 'z_max': None, 't_min': 0, 't_max': 1}
+
+        Rank 1:
+            {'x_min': -91, 'x_max': 50, 'y_min': 475, 'y_max': 500, 'z_min': None, 'z_max': None, 't_min': 0, 't_max': 1}
+
+        Rank 2:
+            {'x_min': -91, 'x_max': 50, 'y_min': 500, 'y_max': 525, 'z_min': None, 'z_max': None, 't_min': 0, 't_max': 1}
+
+        Rank 3:
+            {'x_min': -91, 'x_max': 50, 'y_min': 525, 'y_max': 549, 'z_min': None, 'z_max': None, 't_min': 0, 't_max': 1}
+
+        Calling `get_global_limits` on rank 0 will return:
+
+        >>> {'x_min': -91, 'x_max': 50, 'y_min': 450, 'y_max': 549, 'z_min': None, 'z_max': None, 't_min': 0, 't_max': 1}
+        """
+        my_limits = deepcopy(limits)
+        for key in my_limits:
+            if my_limits[key] is None:
+                my_limits[key] = float('inf') if 'min' in key else float('-inf')
+
+        # Perform MPI reduction to find global min and max values
+        global_limits = {}
+        for key in my_limits:
+            if 'min' in key:
+                global_limits[key] = self.comm.reduce(my_limits[key], op=MPI.MIN, root=0)
+            elif 'max' in key:
+                global_limits[key] = self.comm.reduce(my_limits[key], op=MPI.MAX, root=0)
+
+        # Restore None for dimensions that were originally None
+        if self.rank == 0:
+            for key in global_limits:
+                if global_limits[key] in [float('inf'), float('-inf')]:
+                    global_limits[key] = None
+
+        return global_limits
+
     def _filter_coordinates_selection(self):
         """
         Use the selection limits to filter time, lev, lat, lon, lon_bnds and lat_bnds.
         """
-
-        idx = self._get_idx_intervals()
+        global_limits = self.get_global_limits(self.read_axis_limits)
 
         if self.master:
-            self._full_time = self._full_time[idx["idx_t_min"]:idx["idx_t_max"]]
-            self._full_lev["data"] = self._full_lev["data"][idx["idx_z_min"]:idx["idx_z_max"]]
+            self._full_time = self._full_time[global_limits["t_min"]:global_limits["t_max"]]
+            self._full_lev["data"] = self._full_lev["data"][global_limits["z_min"]:global_limits["z_max"]]
 
             if len(self._full_lat["data"].shape) == 1:
                 # Regular projection
-                self._full_lat["data"] = self._full_lat["data"][idx["idx_y_min"]:idx["idx_y_max"]]
-                self._full_lon["data"] = self._full_lon["data"][idx["idx_x_min"]:idx["idx_x_max"]]
+                self._full_lat["data"] = self._full_lat["data"][global_limits["y_min"]:global_limits["y_max"]]
+                if global_limits["x_min"] < 0:
+                    self._full_lon["data"] = concatenate((self._full_lon["data"][global_limits["x_min"]:],
+                                                          self._full_lon["data"][:global_limits["x_max"]]))
+                    self._full_lon["data"][self._full_lon["data"] > 180] -= 360
+                else:
+                    self._full_lon["data"] = self._full_lon["data"][global_limits["x_min"]:global_limits["x_max"]]
 
                 if self._full_lat_bnds is not None:
-                    self._full_lat_bnds["data"] = self._full_lat_bnds["data"][idx["idx_y_min"]:idx["idx_y_max"], :]
+                    self._full_lat_bnds["data"] = self._full_lat_bnds["data"][global_limits["y_min"]:global_limits["y_max"], :]
                 if self._full_lon_bnds is not None:
-                    self._full_lon_bnds["data"] = self._full_lon_bnds["data"][idx["idx_x_min"]:idx["idx_x_max"], :]
+                    if global_limits["x_min"] < 0:
+                        self._full_lon_bnds["data"] = concatenate((self._full_lon_bnds["data"][global_limits["x_min"]:, :],
+                                                                   self._full_lon_bnds["data"][:global_limits["x_max"], :]))
+
+                        self._full_lon_bnds["data"][self._full_lon_bnds["data"] > 180] -= 360
+                    else:
+                        self._full_lon_bnds["data"] = self._full_lon_bnds["data"][global_limits["x_min"]:global_limits["x_max"], :]
             else:
                 # Irregular projections
-                self._full_lat["data"] = self._full_lat["data"][idx["idx_y_min"]:idx["idx_y_max"],
-                                                                idx["idx_x_min"]:idx["idx_x_max"]]
-                self._full_lon["data"] = self._full_lon["data"][idx["idx_y_min"]:idx["idx_y_max"],
-                                                                idx["idx_x_min"]:idx["idx_x_max"]]
+                if global_limits["x_min"] < 0:
+                    self._full_lat["data"] = concatenate((
+                        self._full_lat["data"][global_limits["y_min"]:global_limits["y_max"],
+                                               global_limits["x_min"]:],
+                        self._full_lat["data"][global_limits["y_min"]:global_limits["y_max"],
+                                               :global_limits["x_max"]]),
+                    axis=1)
+
+                    self._full_lon["data"] = concatenate((
+                        self._full_lon["data"][global_limits["y_min"]:global_limits["y_max"],
+                                               global_limits["x_min"]:],
+                        self._full_lon["data"][global_limits["y_min"]:global_limits["y_max"],
+                                               :global_limits["x_max"]]),
+                    axis=1)
+                    self._full_lon["data"][self._full_lon["data"] > 180] -= 360
+
+                    if self._full_lat_bnds is not None:
+                        self._full_lat_bnds["data"] = concatenate((
+                            self._full_lat_bnds["data"][global_limits["y_min"]:global_limits["y_max"],
+                                                        global_limits["x_min"]:, :],
+                            self._full_lat_bnds["data"][global_limits["y_min"]:global_limits["y_max"],
+                                                        :global_limits["x_max"], :]),
+                    axis=1)
+                    if self._full_lon_bnds is not None:
+                        self._full_lon_bnds["data"] = concatenate((
+                            self._full_lon_bnds["data"][global_limits["y_min"]:global_limits["y_max"],
+                                                        global_limits["x_min"]:, :],
+                            self._full_lon_bnds["data"][global_limits["y_min"]:global_limits["y_max"],
+                                                        :global_limits["x_max"], :]),
+                            axis=1)
+                        self._full_lon_bnds["data"][self._full_lon_bnds["data"] > 180] -= 360
+                else:
+                    self._full_lat["data"] = self._full_lat["data"][global_limits["y_min"]:global_limits["y_max"],
+                                                                    global_limits["x_min"]:global_limits["x_max"]]
 
-                if self._full_lat_bnds is not None:
-                    self._full_lat_bnds["data"] = self._full_lat_bnds["data"][idx["idx_y_min"]:idx["idx_y_max"],
-                                                                              idx["idx_x_min"]:idx["idx_x_max"], :]
-                if self._full_lon_bnds is not None:
-                    self._full_lon_bnds["data"] = self._full_lon_bnds["data"][idx["idx_y_min"]:idx["idx_y_max"],
-                                                                              idx["idx_x_min"]:idx["idx_x_max"], :]
+                    self._full_lon["data"] = self._full_lon["data"][global_limits["y_min"]:global_limits["y_max"],
+                                                                    global_limits["x_min"]:global_limits["x_max"]]
+
+                    if self._full_lat_bnds is not None:
+                        self._full_lat_bnds["data"] = self._full_lat_bnds["data"][global_limits["y_min"]:global_limits["y_max"],
+                                                                                  global_limits["x_min"]:global_limits["x_max"], :]
+                    if self._full_lon_bnds is not None:
+                        self._full_lon_bnds["data"] = self._full_lon_bnds["data"][global_limits["y_min"]:global_limits["y_max"],
+                                                                                  global_limits["x_min"]:global_limits["x_max"], :]
 
         self.hours_start = 0
         self.hours_end = 0
@@ -1507,7 +1624,6 @@ class Nes(object):
             idx["idx_y_max"] = idx_aux
 
         # Axis X
-
         if self.lon_min is None:
             idx["idx_x_min"] = 0
         else:
@@ -1515,7 +1631,11 @@ class Nes(object):
                 axis = 0
             else:
                 axis = 1
-            idx["idx_x_min"] = self._get_coordinate_id(full_lon["data"], self.lon_min, axis=axis)
+            if self.lon_min < 0 and full_lon["data"].max() > 180:
+                aux_xmin = self._get_coordinate_id(full_lon["data"], self.lon_min + 360, axis=axis)
+                idx["idx_x_min"] = - (self.lon["data"].shape[-1] - aux_xmin)
+            else:
+                idx["idx_x_min"] = self._get_coordinate_id(full_lon["data"], self.lon_min, axis=axis)
         if self.lon_max is None:
             idx["idx_x_max"] = full_lon["data"].shape[-1]
         else:
@@ -2454,10 +2574,25 @@ class Nes(object):
                     dim=values["data"].shape))
         elif coordinate_axis == "X":
             if coordinate_len == 1:
-                values["data"] = values["data"][self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"]]
+                if self.read_axis_limits["x_min"] < 0:
+                    # Negative longitudes
+                    values["data"] = concatenate((values["data"][self.read_axis_limits["x_min"]:],
+                                                  values["data"][:self.read_axis_limits["x_max"]]))
+                    values["data"][values["data"] > 180] -= 360
+                else:
+                    values["data"] = values["data"][self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"]]
             elif coordinate_len == 2:
-                values["data"] = values["data"][self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
-                                                self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"]]
+                if self.read_axis_limits["x_min"] < 0:
+                    # Negative longitudes
+                    values["data"] = concatenate(
+                        (values["data"][self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                                               self.read_axis_limits["x_min"]:],
+                                values["data"][self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                                               :self.read_axis_limits["x_max"]]), axis=1)
+                    values["data"][values["data"] > 180] -= 360
+                else:
+                    values["data"] = values["data"][self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                                                    self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"]]
             else:
                 raise NotImplementedError("The coordinate has wrong dimensions: {dim}".format(
                     dim=values["data"].shape))
@@ -2569,55 +2704,123 @@ class Nes(object):
         var_dims = nc_var.dimensions
 
         # Read data in 4 dimensions
-        if len(var_dims) < 2:
-            data = nc_var[:]
-        elif len(var_dims) == 2:
-            data = nc_var[self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
-                          self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"]]
-            data = data.reshape(1, 1, data.shape[-2], data.shape[-1])
-        elif len(var_dims) == 3:
-            if "strlen" in var_dims:
-                data = nc_var[self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
-                              self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"],
-                              :]
-                data_aux = empty(shape=(data.shape[0], data.shape[1]), dtype=object)
-                for lat_n in range(data.shape[0]):
-                    for lon_n in range(data.shape[1]):
-                        data_aux[lat_n, lon_n] = "".join(
-                            data[lat_n, lon_n].tobytes().decode("ascii").replace("\x00", ""))
-                data = data_aux.reshape((1, 1, data_aux.shape[-2], data_aux.shape[-1]))
+        if self.read_axis_limits["x_min"] < 0:
+            if len(var_dims) < 2:
+                data = nc_var[:]
+            elif len(var_dims) == 2:
+                data = concatenate((
+                    nc_var[self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                           self.read_axis_limits["x_min"]:],
+                    nc_var[self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                           :self.read_axis_limits["x_max"]]), axis=1)
+                data = data.reshape(1, 1, data.shape[-2], data.shape[-1])
+            elif len(var_dims) == 3:
+                if "strlen" in var_dims:
+                    data = concatenate((
+                        nc_var[self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                               self.read_axis_limits["x_min"]:, :],
+                        nc_var[self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                               :self.read_axis_limits["x_max"], :]), axis=1)
+                    data_aux = empty(shape=(data.shape[0], data.shape[1]), dtype=object)
+                    for lat_n in range(data.shape[0]):
+                        for lon_n in range(data.shape[1]):
+                            data_aux[lat_n, lon_n] = "".join(
+                                data[lat_n, lon_n].tobytes().decode("ascii").replace("\x00", ""))
+                    data = data_aux.reshape((1, 1, data_aux.shape[-2], data_aux.shape[-1]))
+                else:
+                    data = concatenate((
+                        nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
+                               self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                               self.read_axis_limits["x_min"]:],
+                        nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
+                               self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                               :self.read_axis_limits["x_max"]]), axis=2)
+                    data = data.reshape(data.shape[-3], 1, data.shape[-2], data.shape[-1])
+            elif len(var_dims) == 4:
+                data = concatenate((
+                    nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
+                           self.read_axis_limits["z_min"]:self.read_axis_limits["z_max"],
+                           self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                           self.read_axis_limits["x_min"]:],
+                    nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
+                           self.read_axis_limits["z_min"]:self.read_axis_limits["z_max"],
+                           self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                           :self.read_axis_limits["x_max"]]), axis=3)
+            elif len(var_dims) == 5:
+                if "strlen" in var_dims:
+                    data = concatenate((
+                        nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
+                        self.read_axis_limits["z_min"]:self.read_axis_limits["z_max"],
+                        self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                        self.read_axis_limits["x_min"]:, :],
+                        nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
+                        self.read_axis_limits["z_min"]:self.read_axis_limits["z_max"],
+                        self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                        :self.read_axis_limits["x_max"], :]), axis=3)
+                    data_aux = empty(shape=(data.shape[0], data.shape[1], data.shape[2], data.shape[3]), dtype=object)
+                    for time_n in range(data.shape[0]):
+                        for lev_n in range(data.shape[1]):
+                            for lat_n in range(data.shape[2]):
+                                for lon_n in range(data.shape[3]):
+                                    data_aux[time_n, lev_n, lat_n, lon_n] = "".join(
+                                        data[time_n, lev_n, lat_n, lon_n].tobytes().decode("ascii").replace("\x00", ""))
+                    data = data_aux
+                else:
+                    raise NotImplementedError("Error with {0}. Only can be read netCDF with 4 dimensions or less".format(
+                        var_name))
             else:
-                data = nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
-                              self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                raise NotImplementedError("Error with {0}. Only can be read netCDF with 4 dimensions or less".format(
+                    var_name))
+        else:
+            if len(var_dims) < 2:
+                data = nc_var[:]
+            elif len(var_dims) == 2:
+                data = nc_var[self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
                               self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"]]
-                data = data.reshape(data.shape[-3], 1, data.shape[-2], data.shape[-1])
-        elif len(var_dims) == 4:
-            data = nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
-                          self.read_axis_limits["z_min"]:self.read_axis_limits["z_max"],
-                          self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
-                          self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"]]
-        elif len(var_dims) == 5:
-            if "strlen" in var_dims:
+                data = data.reshape(1, 1, data.shape[-2], data.shape[-1])
+            elif len(var_dims) == 3:
+                if "strlen" in var_dims:
+                    data = nc_var[self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                                  self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"],
+                                  :]
+                    data_aux = empty(shape=(data.shape[0], data.shape[1]), dtype=object)
+                    for lat_n in range(data.shape[0]):
+                        for lon_n in range(data.shape[1]):
+                            data_aux[lat_n, lon_n] = "".join(
+                                data[lat_n, lon_n].tobytes().decode("ascii").replace("\x00", ""))
+                    data = data_aux.reshape((1, 1, data_aux.shape[-2], data_aux.shape[-1]))
+                else:
+                    data = nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
+                                  self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                                  self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"]]
+                    data = data.reshape(data.shape[-3], 1, data.shape[-2], data.shape[-1])
+            elif len(var_dims) == 4:
                 data = nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
                               self.read_axis_limits["z_min"]:self.read_axis_limits["z_max"],
                               self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
-                              self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"],
-                              :]
-                data_aux = empty(shape=(data.shape[0], data.shape[1], data.shape[2], data.shape[3]), dtype=object)
-                for time_n in range(data.shape[0]):
-                    for lev_n in range(data.shape[1]):
-                        for lat_n in range(data.shape[2]):
-                            for lon_n in range(data.shape[3]):
-                                data_aux[time_n, lev_n, lat_n, lon_n] = "".join(
-                                    data[time_n, lev_n, lat_n, lon_n].tobytes().decode("ascii").replace("\x00", ""))
-                data = data_aux
+                              self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"]]
+            elif len(var_dims) == 5:
+                if "strlen" in var_dims:
+                    data = nc_var[self.read_axis_limits["t_min"]:self.read_axis_limits["t_max"],
+                                  self.read_axis_limits["z_min"]:self.read_axis_limits["z_max"],
+                                  self.read_axis_limits["y_min"]:self.read_axis_limits["y_max"],
+                                  self.read_axis_limits["x_min"]:self.read_axis_limits["x_max"],
+                                  :]
+                    data_aux = empty(shape=(data.shape[0], data.shape[1], data.shape[2], data.shape[3]), dtype=object)
+                    for time_n in range(data.shape[0]):
+                        for lev_n in range(data.shape[1]):
+                            for lat_n in range(data.shape[2]):
+                                for lon_n in range(data.shape[3]):
+                                    data_aux[time_n, lev_n, lat_n, lon_n] = "".join(
+                                        data[time_n, lev_n, lat_n, lon_n].tobytes().decode("ascii").replace("\x00", ""))
+                    data = data_aux
+                else:
+                    raise NotImplementedError("Error with {0}. Only can be read netCDF with 4 dimensions or less".format(
+                        var_name))
             else:
                 raise NotImplementedError("Error with {0}. Only can be read netCDF with 4 dimensions or less".format(
                     var_name))
-        else:
-            raise NotImplementedError("Error with {0}. Only can be read netCDF with 4 dimensions or less".format(
-                var_name))
-
+        
         # Unmask array
         data = self._unmask_array(data)
 
diff --git a/setup.py b/setup.py
index 76bc9dc5529c5a07f4028c18238b0f8ed7d64188..0e300345b547a09a2330b8f0997fd04a1003e53a 100755
--- a/setup.py
+++ b/setup.py
@@ -34,7 +34,7 @@ setup(
     name='nes',
     license='Apache License 2.0',
     version=__version__,
-    description='',
+    description='NES (NetCDF for Earth Science) is a library designed to operate, read, and write NetCDF datasets in a parallel and efficient yet transparent way for the user, using domain decomposition as the parallelization strategy with customizable axis splitting.',
     long_description=long_description,
     long_description_content_type="text/markdown",
     author="Carles Tena Medina, Alba Vilanova Cortezón, Carmen Piñero Megías",
diff --git a/tests/1.4-test_selecting_negatives.py b/tests/1.4-test_selecting_negatives.py
new file mode 100644
index 0000000000000000000000000000000000000000..2c0e4ec995b97480445225d391d0d74e64af8f46
--- /dev/null
+++ b/tests/1.4-test_selecting_negatives.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+import sys
+import timeit
+import pandas as pd
+from mpi4py import MPI
+from nes import open_netcdf
+from datetime import datetime
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+size = comm.Get_size()
+
+parallel_method = 'Y'
+serial_write = True
+
+result_path = "Times_test_1.4.Selecting_{0}_{1:03d}.csv".format(parallel_method, size)
+
+result = pd.DataFrame(index=['read', 'calcul', 'write'],
+                      columns=['1.4.1.Selecting_Negatives'])
+
+# NAMEE
+src_path = "/gpfs/projects/bsc32/models/NES_tutorial_data/ga_20220101.nc"
+var_list = ['so2', "pm25", 'nox']
+
+# ======================================================================================================================
+# ======================================      '1.4.1.Selecting_Negatives'       =====================================================
+# ======================================================================================================================
+test_name = '1.4.1.Selecting_Negatives'
+
+if rank == 0:
+    print(test_name)
+
+st_time = timeit.default_timer()
+
+# Source data
+nessy = open_netcdf(src_path, parallel_method=parallel_method, balanced=True)
+nessy.keep_vars(var_list)
+nessy.sel(lat_min=35, lat_max=45, lon_min=-9, lon_max=5)
+nessy.cell_measures = {}
+# nessy.calculate_grid_area(overwrite=True)
+nessy.load()
+print(f"Rank {nessy.rank} -> Lon: {nessy.lon}")
+print(f"Rank {nessy.rank} -> Lon BNDS: {nessy.lon_bnds}")
+print(f"Rank {nessy.rank} -> FULL Lon: {nessy.get_full_longitudes()}")
+
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+nessy.to_netcdf(test_name.replace(' ', '_') + "{0:03d}.nc".format(size), serial=serial_write)
+
+comm.Barrier()
+result.loc['write', test_name] = timeit.default_timer() - st_time
+
+comm.Barrier()
+if rank == 0:
+    print(result.loc[:, test_name])
+sys.stdout.flush()
+
+if rank == 0:
+    result.to_csv(result_path)
+    print("TEST PASSED SUCCESSFULLY!!!!!")
diff --git a/tests/test_bash.mn5.sh b/tests/test_bash.mn5.sh
index e767c3e402ff091e317de78211364db49a2aba03..606abec246c19d8866c34933a30d84fd636349b7 100644
--- a/tests/test_bash.mn5.sh
+++ b/tests/test_bash.mn5.sh
@@ -15,11 +15,14 @@ source activate base
 conda activate /gpfs/projects/bsc32/repository/apps/conda_envs/NES_v1.1.9
 
 
+
 cd /gpfs/scratch/bsc32/${USER}/AC_PostProcess/NES/tests || exit
 
 mpirun -np 4 python 1.1-test_read_write_projection.py
 mpirun -np 4 python 1.2-test_create_projection.py
 mpirun -np 4 python 1.3-test_selecting.py
+mpirun -np 4 python 1.4-test_selecting_negatives.py
+
 
 mpirun -np 4 python 2.1-test_spatial_join.py
 mpirun -np 4 python 2.2-test_create_shapefile.py