diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py
index 7b0759831cd17644ac85f83e64e03354738acea3..2d056372d0a2a88377f77bd0be6b94b870d3ef24 100644
--- a/nes/nc_projections/default_nes.py
+++ b/nes/nc_projections/default_nes.py
@@ -9,7 +9,6 @@ from xarray import open_dataset
 from netCDF4 import Dataset, num2date, date2num
 from mpi4py import MPI
 from cfunits import Units
-from numpy.ma.core import MaskError
 import geopandas as gpd
 from shapely.geometry import Polygon, Point
 from copy import deepcopy, copy
@@ -69,7 +68,7 @@ class Nes(object):
     write_axis_limits : dict
         Dictionary with the 4D limits of the rank data to write.
         t_min, t_max, z_min, z_max, y_min, y_max, x_min and x_max.
-    time : List
+    time : List[datetime]
         List of time steps of the rank data.
     lev : dict
         Vertical levels dictionary with the portion of 'data' corresponding to the rank values.
@@ -79,12 +78,12 @@ class Nes(object):
         Longitudes dictionary with the portion of 'data' corresponding to the rank values.
     global_attrs : dict
         Global attributes with the attribute name as key and data as values.
-    _var_dim : None, tuple
-        Tuple with the name of the Y and X dimensions for the variables.
-    _lat_dim : None, tuple
-        Tuple with the name of the dimensions of the Latitude values.
-    _lon_dim : None, tuple
-        Tuple with the name of the dimensions of the Longitude values.
+    _var_dim : None or tuple
+        Name of the Y and X dimensions for the variables.
+    _lat_dim : None or tuple
+        Name of the dimensions of the Latitude values.
+    _lon_dim : None or tuple
+        Name of the dimensions of the Longitude values.
     """
     def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
@@ -102,7 +101,7 @@ class Nes(object):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
+        xarray: bool
             (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default over Y axis
@@ -116,11 +115,11 @@ class Nes(object):
             Number of hours to remove from last time steps.
         first_level : int
             Index of the first level to use.
-        last_level : int, None
+        last_level : int or None
             Index of the last level to use. None if it is the last.
         create_nes : bool
             Indicates if you want to create the object from scratch (True) or through an existing file.
-        times : List, None
+        times : List[datetime] or None
             List of times to substitute the current ones while creation.
         kwargs :
             Projection dependent parameters to create it from scratch
@@ -272,7 +271,7 @@ class Nes(object):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
+        xarray: bool
             (Not working) Indicates if you want to use xarray as default.
         avoid_first_hours : int
             Number of hours to remove from first time steps.
@@ -286,7 +285,7 @@ class Nes(object):
             Balanced dataset cannot be written in chunking mode.
         first_level : int
             Index of the first level to use.
-        last_level : int, None
+        last_level : int or None
             Index of the last level to use. None if it is the last.
         create_nes : bool
             Indicates if you want to create the object from scratch (True) or through an existing file.
@@ -364,6 +363,35 @@ class Nes(object):
 
         return None
 
+    def __add__(self, other):
+        """
+        Sum two NES objects
+
+        Parameters
+        ----------
+        other : Nes
+            Nes to be summed
+
+        Returns
+        -------
+        Nes
+            Summed Nes object
+        """
+        nessy = self.copy(copy_vars=True)
+        for var_name in other.variables.keys():
+            if var_name not in nessy.variables.keys():
+                # Create New variable
+                nessy.variables[var_name] = deepcopy(other.variables[var_name])
+            else:
+                nessy.variables[var_name]['data'] += other.variables[var_name]['data']
+        return nessy
+
+    def __radd__(self, other):
+        if other == 0 or other is None:
+            return self
+        else:
+            return self.__add__(other)
+
     def copy(self, copy_vars=False):
         """
         Copy the Nes object.
@@ -383,7 +411,8 @@ class Nes(object):
         nessy = deepcopy(self)
         nessy.netcdf = None
         if copy_vars:
-            nessy.variables = nessy._get_lazy_variables()
+            nessy.set_communicator(self.comm)
+            nessy.variables = deepcopy(self.variables)
             nessy.cell_measures = deepcopy(self.cell_measures)
         else:
             nessy.variables = {}
@@ -494,7 +523,8 @@ class Nes(object):
 
         return None
 
-    def create_single_spatial_bounds(self, coordinates, inc, spatial_nv=2, inverse=False):
+    @staticmethod
+    def create_single_spatial_bounds(coordinates, inc, spatial_nv=2, inverse=False):
         """
         Calculate the vertices coordinates.
 
@@ -505,7 +535,7 @@ class Nes(object):
         inc : float 
             Increment between centre values.
         spatial_nv : int 
-            Non mandatory parameter that informs the number of vertices that the boundaries must have. Default: 2.
+            Non-mandatory parameter that informs the number of vertices that the boundaries must have. Default: 2.
         inverse : bool 
             For some grid latitudes.
 
@@ -544,20 +574,14 @@ class Nes(object):
         inc_lat = np.abs(np.mean(np.diff(self._lat['data'])))
         lat_bnds = self.create_single_spatial_bounds(self._lat['data'], inc_lat, spatial_nv=2)
 
-        self._lat_bnds = {}
-        self._lat_bnds['data'] = deepcopy(lat_bnds)
-        self.lat_bnds = {}
-        self.lat_bnds['data'] = lat_bnds[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'],
-                                         :]
+        self._lat_bnds = {'data': deepcopy(lat_bnds)}
+        self.lat_bnds = {'data': lat_bnds[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'], :]}
         
         inc_lon = np.abs(np.mean(np.diff(self._lon['data'])))
         lon_bnds = self.create_single_spatial_bounds(self._lon['data'], inc_lon, spatial_nv=2)
 
-        self._lon_bnds = {}
-        self._lon_bnds['data'] = deepcopy(lon_bnds)
-        self.lon_bnds = {}
-        self.lon_bnds['data'] = lon_bnds[self.read_axis_limits['x_min']:self.read_axis_limits['x_max'], 
-                                         :]
+        self._lon_bnds = {'data': deepcopy(lon_bnds)}
+        self.lon_bnds = {'data': lon_bnds[self.read_axis_limits['x_min']:self.read_axis_limits['x_max'], :]}
 
         return None
 
@@ -569,9 +593,9 @@ class Nes(object):
 
         Returns
         -------
-        lon_bnds_mesh : numpy.array
+        lon_bnds_mesh : numpy.ndarray
             Longitude boundaries in the mesh format
-        lat_bnds_mesh : numpy.array
+        lat_bnds_mesh : numpy.ndarray
             Latitude boundaries in the mesh format
         """
         if self.size > 1:
@@ -601,17 +625,17 @@ class Nes(object):
             lon_bnds_mesh[1:, 1:] = self.lon_bnds['data'][:, :, 2]
             lon_bnds_mesh[1:, :-1] = self.lon_bnds['data'][:, :, 3]
         else:
-            raise RuntimeError("Invalid number of vertices: {0}".format(self.lat_bnds['data'].shape[-1] ))
+            raise RuntimeError("Invalid number of vertices: {0}".format(self.lat_bnds['data'].shape[-1]))
 
         return lon_bnds_mesh, lat_bnds_mesh
 
     def free_vars(self, var_list):
         """
-        Erase the selected variables from the variables information.
+        Erase the selected variables from the variables' information.
 
         Parameters
         ----------
-        var_list : List, str, list
+        var_list : List or str
             List (or single string) of the variables to be loaded.
         """
 
@@ -638,7 +662,7 @@ class Nes(object):
 
         Parameters
         ----------
-        var_list : List, str
+        var_list : List or str
             List (or single string) of the variables to be loaded.
         """
 
@@ -666,7 +690,7 @@ class Nes(object):
 
         return time_interval
 
-    def sel_time(self, time, copy=False):
+    def sel_time(self, time, do_copy=False):
         """
         To select only one time step.
 
@@ -674,7 +698,7 @@ class Nes(object):
         ----------
         time : datetime.datetime
             Time stamp to select.
-        copy : bool
+        do_copy : bool
             Indicates if you want a copy with the selected time step (True) or to modify te existing one (False).
 
         Returns
@@ -683,7 +707,7 @@ class Nes(object):
             Nes object with the data (and metadata) of the selected time step.
         """
 
-        if copy:
+        if do_copy:
             aux_nessy = self.copy(copy_vars=False)
             aux_nessy.comm = self.comm
         else:
@@ -982,7 +1006,7 @@ class Nes(object):
             self.set_time_bnds(aux_time_bounds)
 
         elif type_op == 'withoutt0':
-            for var_name, var_info in self.variables.items ():
+            for var_name, var_info in self.variables.items():
                 if var_info['data'] is None:
                     self.load(var_name)
                 if op == 'mean':
@@ -1197,10 +1221,10 @@ class Nes(object):
 
         rows_sum = 0
         for proc in range(self.size):
-            fid_dist[proc] = {'x_min': None, 'x_max': None,
-                              'y_min': None, 'y_max': None,
-                              'z_min': None, 'z_max': None,
-                              't_min': None, 't_max': None}
+            fid_dist[proc] = {'x_min': 0, 'x_max': None,
+                              'y_min': 0, 'y_max': None,
+                              'z_min': 0, 'z_max': None,
+                              't_min': 0, 't_max': None}
             if proc < procs_rows_extended:
                 aux_rows = procs_len + 1
             else:
@@ -1542,9 +1566,9 @@ class Nes(object):
         Returns
         -------
         lat_bnds : List
-            List of latitude bounds of the NetCDF data.
+            Latitude bounds of the NetCDF data.
         lon_bnds : List
-            List of longitude bounds of the NetCDF data.
+            Longitude bounds of the NetCDF data.
         """
 
         if self.is_xarray:
@@ -1554,13 +1578,11 @@ class Nes(object):
             if self.master:
                 if not create_nes:
                     if 'lat_bnds' in self.netcdf.variables.keys():
-                        lat_bnds = {}
-                        lat_bnds['data'] = self.netcdf.variables['lat_bnds'][:]
+                        lat_bnds = {'data': self.netcdf.variables['lat_bnds'][:]}
                     else:
                         lat_bnds = None
                     if 'lon_bnds' in self.netcdf.variables.keys():
-                        lon_bnds = {}
-                        lon_bnds['data'] = self.netcdf.variables['lon_bnds'][:]
+                        lon_bnds = {'data': self.netcdf.variables['lon_bnds'][:]}
                     else:
                         lon_bnds = None
                 else:
@@ -1587,21 +1609,21 @@ class Nes(object):
 
         Returns
         -------
-        cell_measures : dict
+        dict
             Dictionary of cell measures of the NetCDF data.
         """
 
-        cell_measures = {}
+        c_measures = {}
         if self.master:
             if not create_nes:
                 if 'cell_area' in self.netcdf.variables.keys():
-                    cell_measures['cell_area'] = {}
-                    cell_measures['cell_area']['data'] = self.netcdf.variables['cell_area'][:]
-        cell_measures = self.comm.bcast(cell_measures, root=0)
+                    c_measures['cell_area'] = {}
+                    c_measures['cell_area']['data'] = self.netcdf.variables['cell_area'][:]
+        c_measures = self.comm.bcast(c_measures, root=0)
 
         self.free_vars(['cell_area'])
 
-        return cell_measures
+        return c_measures
 
     def _get_coordinate_dimension(self, possible_names):
         """
@@ -1789,6 +1811,7 @@ class Nes(object):
         data: np.array
             Portion of the variable data corresponding to the rank.
         """
+        # from numpy.ma.core import MaskError
 
         nc_var = self.netcdf.variables[var_name]
         var_dims = nc_var.dimensions
@@ -1810,7 +1833,7 @@ class Nes(object):
                     for lon_n in range(data.shape[1]):
                         data_aux[lat_n, lon_n] = ''.join(
                             data[lat_n, lon_n].tostring().decode('ascii').replace('\x00', ''))
-                data = data_aux.reshape(1, 1, data_aux.shape[-2], data_aux.shape[-1])
+                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'],
@@ -1830,11 +1853,11 @@ class Nes(object):
         else:
             raise NotImplementedError('Error with {0}. Only can be read netCDF with 4 dimensions or less'.format(
                 var_name))
-        # Missing to nan
-        try:
-            data[data.shapefile == True] = np.nan
-        except (AttributeError, MaskError, ValueError):
-            pass
+        # # Missing to nan
+        # try:
+        #     data[data.shapefile == True] = np.nan
+        # except (AttributeError, MaskError, ValueError):
+        #     pass
 
         return data
 
@@ -2042,10 +2065,10 @@ class Nes(object):
 
         rows_sum = 0
         for proc in range(self.size):
-            fid_dist[proc] = {'x_min': None, 'x_max': None,
-                              'y_min': None, 'y_max': None,
-                              'z_min': None, 'z_max': None,
-                              't_min': None, 't_max': None}
+            fid_dist[proc] = {'x_min': 0, 'x_max': None,
+                              'y_min': 0, 'y_max': None,
+                              'z_min': 0, 'z_max': None,
+                              't_min': 0, 't_max': None}
             if proc < procs_rows_extended:
                 aux_rows = procs_len + 1
             else:
@@ -2109,7 +2132,7 @@ class Nes(object):
 
         # TIMES
         time_var = netcdf.createVariable('time', np.float64, ('time',), zlib=self.zip_lvl > 0, complevel=self.zip_lvl)
-        time_var.units = 'hours since {0}'.format( self._time[0].strftime('%Y-%m-%d %H:%M:%S'))
+        time_var.units = 'hours since {0}'.format(self._time[0].strftime('%Y-%m-%d %H:%M:%S'))
         time_var.standard_name = 'time'
         time_var.calendar = 'standard'
         time_var.long_name = 'time'
@@ -2217,7 +2240,7 @@ class Nes(object):
             if var_dict['data'] is not None:
                 
                 # Get dimensions
-                if var_dict['data'].shape == 4:
+                if len(var_dict['data'].shape) == 4:
                     var_dims = ('time', 'lev',) + self._var_dim
                 else:
                     var_dims = self._var_dim
@@ -2228,13 +2251,13 @@ class Nes(object):
                     if var_dtype != var_dict['data'].dtype:
                         msg = "WARNING!!! "
                         msg += "Different data types for variable {0}. ".format(var_name)
-                        msg += "Input dtype={0}. Data dtype={1}.".format(var_dtype, 
-                                                                        var_dict['data'].dtype)
+                        msg += "Input dtype={0}. Data dtype={1}.".format(var_dtype, var_dict['data'].dtype)
                         warnings.warn(msg)
                         try:
                             var_dict['data'] = var_dict['data'].astype(var_dtype)
                         except Exception as e:  # TODO: Detect exception
-                            raise e("It was not possible to cast the data to the input dtype.")
+                            print(e)
+                            raise TypeError("It was not possible to cast the data to the input dtype.")
                 else:
                     var_dtype = var_dict['data'].dtype
 
@@ -2244,27 +2267,37 @@ class Nes(object):
                     var_dtype = var_dict['data'].dtype
 
                 # Convert list of strings to chars for parallelization
-                try:
-                    unicode_type = len(max(var_dict['data'].flatten(), key=len))
-                    if ((var_dict['data'].dtype == np.dtype('<U' + str(unicode_type))) or 
-                        (var_dict['data'].dtype == np.str)):
-                        # Add strlen as a dimension if needed
-                        if 'strlen' not in var_dims:
-                            var_dims += ('strlen',)
-                        data_aux = np.chararray(shape=(var_dict['data'].shape[0], 
-                                                       var_dict['data'].shape[1], 
-                                                       self.strlen))
-                        for lat_n in range(var_dict['data'].shape[0]):
-                            for lon_n in range(var_dict['data'].shape[1]):
-                                chars = [v for v in var_dict['data'][lat_n][lon_n]]
-                                data = chars + ['']*(self.strlen-len(chars))
-                                for char_n, char in enumerate(data):
-                                    data_aux[lat_n, lon_n, char_n] = char.encode('ascii', 'ignore')
+                if (var_dict['data'].dtype != float) and (var_dict['data'].dtype != int):
+                    try:
+                        # Get unicode 
+                        unicode_type = len(max(var_dict['data'].flatten(), key=len))
                         
-                        var_dict['data'] = data_aux
-                        var_dtype = 'S1'
-                except:
-                    pass
+                        if ((var_dict['data'].dtype == np.dtype('<U' + str(unicode_type))) or
+                            (var_dict['data'].dtype == str)):
+                            
+                            # Add strlen as a dimension if needed
+                            if 'strlen' not in var_dims:
+                                var_dims += ('strlen',)
+
+                            # Split strings into chars (S1)
+                            data_aux = np.chararray(shape=(var_dict['data'].shape[0], 
+                                                           var_dict['data'].shape[1], 
+                                                           self.strlen))
+                            for lat_n in range(var_dict['data'].shape[0]):
+                                for lon_n in range(var_dict['data'].shape[1]):
+                                    chars = [v for v in var_dict['data'][lat_n][lon_n]]
+                                    data = chars + ['']*(self.strlen-len(chars))
+                                    for char_n, char in enumerate(data):
+                                        data_aux[lat_n, lon_n, char_n] = char.encode('ascii', 'ignore')      
+                            var_dict['data'] = data_aux
+                            var_dtype = 'S1'
+
+                    # TODO: Detect exception
+                    except Exception as e:
+                        print('String data of {0} could not be converted into chars while writing'.format(var_name))
+                        print('Error:' + e)
+                        # TODO document that exception
+                        pass
 
                 if self.info:
                     print("Rank {0:03d}: Writing {1} var ({2}/{3})".format(
@@ -2303,14 +2336,15 @@ class Nes(object):
                                         self.write_axis_limits['z_min']:self.write_axis_limits['z_max'],
                                         self.write_axis_limits['y_min']:self.write_axis_limits['y_max'],
                                         self.write_axis_limits['x_min']:self.write_axis_limits['x_max']] = att_value
-                                except ValueError:
-                                    var[self.write_axis_limits['t_min']:self.write_axis_limits['t_max'],
-                                        0,
-                                        self.write_axis_limits['y_min']:self.write_axis_limits['y_max'],
-                                        self.write_axis_limits['x_min']:self.write_axis_limits['x_max']] = att_value
-                                    # msg = "*WARNING* '{0}' variable is a 3D field. Setting it on first (0) layer.".format(
-                                    #     var_name)
-                                    # warn(msg)
+                                except ValueError as e:
+                                    print(var)
+                                    print(att_value)
+                                    raise e
+                                    # var[self.write_axis_limits['t_min']:self.write_axis_limits['t_max'],
+                                    #     0,
+                                    #     self.write_axis_limits['y_min']:self.write_axis_limits['y_max'],
+                                    #     self.write_axis_limits['x_min']:self.write_axis_limits['x_max']] = att_value
+
                                 except IndexError:
                                     raise IndexError("Different shapes. out_shape={0}, data_shp={1}".format(
                                         var[self.write_axis_limits['t_min']:self.write_axis_limits['t_max'],
@@ -2333,15 +2367,15 @@ class Nes(object):
                                 else:
                                     raise NotImplementedError('It is not possible to write 3D variables.')
                             if self.info:
-                                print("Rank {0:03d}: Var {1} data ({2}/{3})".format(self.rank, var_name, i + 1,
-                                                                                    len(self.variables)))
+                                print("Rank {0:03d}: Var {1} data ({2}/{3})".format(
+                                    self.rank, var_name, i + 1, len(self.variables)))
                         elif att_name not in ['chunk_size', 'var_dims', 'dimensions']:
                             var.setncattr(att_name, att_value)
 
                     self._set_var_crs(var)
                     if self.info:
-                        print("Rank {0:03d}: Var {1} completed ({2}/{3})".format(self.rank, var_name, i + 1,
-                                                                                len(self.variables)))
+                        print("Rank {0:03d}: Var {1} completed ({2}/{3})".format(
+                            self.rank, var_name, i + 1, len(self.variables)))
                 except Exception as e:
                     print("**ERROR** an error has occurred while writing the '{0}' variable".format(var_name))
                     # print("**ERROR** an error has occurredred while writing the '{0}' variable".format(var_name),
@@ -2354,28 +2388,25 @@ class Nes(object):
 
         return None
 
-    def _create_centre_coordinates(self):
-        """
-        Must be implemented on inner class.
+    def _create_centre_coordinates(self, **kwargs):
         """
+        Calculate centre latitudes and longitudes from grid details.
 
-        return None
+        Must be implemented on inner classes
 
-    def _create_metadata(self, netcdf):
-        """
-        Must be implemented on inner class.
+        Returns
+        ----------
+        centre_lat : dict
+            Dictionary with data of centre latitudes in 1D
+        centre_lon : dict
+            Dictionary with data of centre longitudes in 1D
         """
 
         return None
 
-    def _set_crs(self, netcdf):
+    def _create_metadata(self, netcdf):
         """
         Must be implemented on inner class.
-
-        Parameters
-        ----------
-        netcdf : Dataset
-            netcdf4-python Dataset.
         """
 
         return None
@@ -2463,9 +2494,11 @@ class Nes(object):
         info : bool
             Indicates if you want to print the information of each writing step by stdout Default: False.
         chunking : bool
-            Indicates if you want a chunked netCDF output. Only available with non serial writes. Default: False.
+            Indicates if you want a chunked netCDF output. Only available with non-serial writes. Default: False.
+        type : str
+            Type to NetCDf to write. 'CAMS_RA' or 'NES'
         """
-
+        nc_type = type
         old_info = self.info
         self.info = info
 
@@ -2488,14 +2521,14 @@ class Nes(object):
                         new_nc.__to_netcdf_cams_ra(path)
                     else:
                         raise ValueError(
-                            "Unknown NetCDF type '{0}'. Use 'CAMS_RA' or 'NES'; default='NES'".format(type))
+                            "Unknown NetCDF type '{0}'. Use 'CAMS_RA' or 'NES'; default='NES'".format(nc_type))
             else:
-                if type == 'NES':
+                if nc_type == 'NES':
                     self.__to_netcdf_py(path, chunking=chunking)
-                elif type == 'CAMS_RA':
+                elif nc_type == 'CAMS_RA':
                     self.__to_netcdf_cams_ra(path)
                 else:
-                    raise ValueError("Unknown NetCDF type '{0}'. Use 'CAMS_RA' or 'NES'; default='NES'".format(type))
+                    raise ValueError("Unknown NetCDF type '{0}'. Use 'CAMS_RA' or 'NES'; default='NES'".format(nc_type))
 
         self.info = old_info
 
@@ -2611,17 +2644,21 @@ class Nes(object):
             Dictionary with the grib2 keys.
         grib_template_path : str
             Path to the grib2 file to use as template.
+        lat_flip : bool
+            Indicates if the latitude values (and data) has to be flipped
         info : bool
             Indicates if you want to print extra information during the process.
         """
 
         # if serial:
         if self.parallel_method in ['X', 'Y'] and self.size > 1:
-            data = self._gather_data()
+            data = self._gather_data(self.variables)
+            c_measures = self._gather_data(self.cell_measures)
             if self.master:
                 new_nc = self.copy(copy_vars=False)
                 new_nc.set_communicator(MPI.COMM_SELF)
                 new_nc.variables = data
+                new_nc.cell_measures = c_measures
                 new_nc.__to_grib2(path, grib_keys, grib_template_path, lat_flip=lat_flip, info=info)
         else:
             self.__to_grib2(path, grib_keys, grib_template_path, lat_flip=lat_flip, info=info)
@@ -2779,8 +2816,8 @@ class Nes(object):
         """
         Add variables data to shapefile.
 
-        var_list : List
-            List (or single string) of the variables to be loaded and saved in the shapefile.
+        var_list : List or str
+           Variables to be loaded and saved in the shapefile.
         idx_lev : int
             Index of vertical level for which the data will be saved in the shapefile.
         idx_time : int
@@ -2798,12 +2835,14 @@ class Nes(object):
 
         Parameters
         ----------
-        ext_shp : GeoPandasDataFrame, str
+        ext_shp : GeoPandasDataFrame or str
             File or path from where the data will be obtained on the intersection.
         method : str
             Overlay method. Accepted values: ['nearest', 'intersection', 'centroid'].
-        var_list : List, None
+        var_list : List or None
             Variables that will be included in the resulting shapefile.
+        info : bool
+            Indicates if you want to print the process info or no
         """
         if isinstance(ext_shp, str):
             ext_shp = gpd.read_file(ext_shp)
@@ -2837,7 +2876,8 @@ class Nes(object):
 
             # From geodetic coordinates (e.g. 4326) to meters (e.g. 4328) to use sjoin_nearest
             # TODO: Check if the projection 4328 does not distort the coordinates too much
-            # https://gis.stackexchange.com/questions/372564/userwarning-when-trying-to-get-centroid-from-a-polygon-geopandas
+            # https://gis.stackexchange.com/questions/372564/
+            #   userwarning-when-trying-to-get-centroid-from-a-polygon-geopandas
             ext_shp = ext_shp.to_crs('EPSG:4328')
             centroids_gdf = centroids_gdf.to_crs('EPSG:4328')
             
@@ -2875,7 +2915,7 @@ class Nes(object):
 
                 intersection.loc[:, 'weight'] = 1.
                 for i, row in intersection.iterrows():
-                    if i % 1000 == 0 and self.master and info:
+                    if isinstance(i, int) and i % 1000 == 0 and self.master and info:
                         print('\tRank {0:03d}: {1:.3f} %'.format(self.rank, i*100 / len(intersection)))
                     # Filter to do not calculate percentages over 100% grid cells spatial joint
                     if counts[row['FID']] > 1:
@@ -2943,7 +2983,7 @@ class Nes(object):
 
         return centroids_gdf
 
-    def __gather_data_py_object(self):
+    def __gather_data_py_object(self, data_to_gather):
         """
         Gather all the variable data into the MPI rank 0 to perform a serial write.
 
@@ -2953,7 +2993,7 @@ class Nes(object):
             Variables dictionary with all the data from all the ranks.
         """
 
-        data_list = deepcopy(self.variables)
+        data_list = deepcopy(data_to_gather)
         for var_name in data_list.keys():
             try:
                 # noinspection PyArgumentList
@@ -3029,12 +3069,13 @@ class Nes(object):
                 print("Gathering {0}".format(var_name))
             shp_len = len(data_list[var_name]['data'].shape)
             try:
-                # Collect local array sizes using the high-level mpi4py gather
+                # Collect local array sizes using the gather communication pattern
                 rank_shapes = np.array(self.comm.gather(data_list[var_name]['data'].shape, root=0))
                 sendbuf = data_list[var_name]['data'].flatten()
                 sendcounts = np.array(self.comm.gather(len(sendbuf), root=0))
                 if self.master:
-                    recvbuf = np.empty(sum(sendcounts), dtype=type(sendbuf[0]))
+                    # recvbuf = np.empty(sum(sendcounts), dtype=type(sendbuf[0]))
+                    recvbuf = np.empty(sum(sendcounts), dtype=type(sendbuf.max()))
                 else:
                     recvbuf = None
                 self.comm.Gatherv(sendbuf=sendbuf, recvbuf=(recvbuf, sendcounts), root=0)
@@ -3151,7 +3192,7 @@ class Nes(object):
         self : Nes
             Source Nes object.
         new_levels : List
-            List of new vertical levels.
+            New vertical levels.
         new_src_vertical
         kind : str
             Vertical methods type.
diff --git a/nes/nc_projections/points_nes.py b/nes/nc_projections/points_nes.py
index 5a2edaa418dae585cdfaaffc6ab6a38f8cb848f3..318302e23519ac3ce5ae0fb9e58353fc8a0b0de2 100644
--- a/nes/nc_projections/points_nes.py
+++ b/nes/nc_projections/points_nes.py
@@ -368,20 +368,32 @@ class PointsNes(Nes):
                             var_dims = ('time',) + self._var_dim
 
                     # Convert list of strings to chars for parallelization
-                    try:
-                        unicode_type = len(max(var_dict['data'].flatten(), key=len))
-                        if ((var_dict['data'].dtype == np.dtype('<U' + str(unicode_type))) or 
-                            (var_dict['data'].dtype == np.str)):
-                            # Add strlen as a dimension if needed
-                            if self.strlen is not None:
-                                if 'strlen' not in var_dims:
-                                    var_dims += ('strlen',)
-                                var_dict['data'] = stringtochar(np.array([v.encode('ascii', 'ignore') 
-                                                                for v in var_dict['data']]).astype('S' + 
-                                                                                                   str(self.strlen)))
-                                var_dtype = 'S1'
-                    except:
-                        pass
+                    if (var_dict['data'].dtype != float) and (var_dict['data'].dtype != int):
+                        try:
+                            # Get unicode 
+                            unicode_type = len(max(var_dict['data'].flatten(), key=len))
+                            
+                            if ((var_dict['data'].dtype == np.dtype('<U' + str(unicode_type))) or 
+                                (var_dict['data'].dtype == np.str)):
+
+                                if self.strlen is not None:
+
+                                    # Add strlen as a dimension if needed
+                                    if 'strlen' not in var_dims:
+                                        var_dims += ('strlen',)
+
+                                    # Split strings into chars (S1)
+                                    var_dict['data'] = stringtochar(np.array([v.encode('ascii', 'ignore') 
+                                                                    for v in var_dict['data']]).astype('S' + 
+                                                                                                    str(self.strlen)))
+                                    var_dtype = 'S1'
+                        
+                        # TODO: Detect exception
+                        except Exception as e:
+                            print('String data of {0} could not be converted into chars while writing'.format(var_name))
+                            print('Error:' + e)
+                            # TODO document that exception
+                            pass
                     
                     # Ensure data is of type numpy array (to create NES)
                     if not isinstance(var_dict['data'], (np.ndarray, np.generic)):
diff --git a/nes/nc_projections/points_nes_ghost.py b/nes/nc_projections/points_nes_ghost.py
index d7cf59891e6355341d3c7e8506dd25e60bf8ba28..695ac4b658dc36b38f9434f980b470eafc2e2ff6 100644
--- a/nes/nc_projections/points_nes_ghost.py
+++ b/nes/nc_projections/points_nes_ghost.py
@@ -350,20 +350,32 @@ class PointsNesGHOST(PointsNes):
                         var_dims = self._var_dim + ('time',)
 
                     # Convert list of strings to chars for parallelization
-                    try:
-                        unicode_type = len(max(var_dict['data'].flatten(), key=len))
-                        if ((var_dict['data'].dtype == np.dtype('<U' + str(unicode_type))) or 
-                            (var_dict['data'].dtype == np.str)):
-                            # Add strlen as a dimension if needed
-                            if self.strlen is not None:
-                                if 'strlen' not in var_dims:
-                                    var_dims += ('strlen',)
-                                var_dict['data'] = stringtochar(np.array([v.encode('ascii', 'ignore') 
-                                                                for v in var_dict['data']]).astype('S' + 
-                                                                                                   str(self.strlen)))
-                                var_dtype = 'S1'
-                    except:
-                        pass
+                    if (var_dict['data'].dtype != float) and (var_dict['data'].dtype != int):
+                        try:
+                            # Get unicode
+                            unicode_type = len(max(var_dict['data'].flatten(), key=len))
+                            
+                            if ((var_dict['data'].dtype == np.dtype('<U' + str(unicode_type))) or 
+                                (var_dict['data'].dtype == np.str)):
+
+                                if self.strlen is not None:
+
+                                    # Add strlen as a dimension if needed
+                                    if 'strlen' not in var_dims:
+                                        var_dims += ('strlen',)
+                                    
+                                    # Split strings into chars (S1)
+                                    var_dict['data'] = stringtochar(np.array([v.encode('ascii', 'ignore') 
+                                                                    for v in var_dict['data']]).astype('S' + 
+                                                                                                    str(self.strlen)))
+                                    var_dtype = 'S1'
+
+                        # TODO: Detect exception
+                        except Exception as e:
+                            print('String data of {0} could not be converted into chars while writing'.format(var_name))
+                            print('Error:' + e)
+                            # TODO document that exception
+                            pass
                     
                     if self.info:
                         print("Rank {0:03d}: Writing {1} var ({2}/{3})".format(self.rank, var_name, i + 1,
diff --git a/nes/nc_projections/points_nes_providentia.py b/nes/nc_projections/points_nes_providentia.py
index 9d576ad2923dad041c7f448513913855003f2fc7..3e45308b688876d935682f1cfff3ac8da66e02b0 100644
--- a/nes/nc_projections/points_nes_providentia.py
+++ b/nes/nc_projections/points_nes_providentia.py
@@ -386,20 +386,32 @@ class PointsNesProvidentia(PointsNes):
                         var_dims = self._var_dim + ('time',)
 
                     # Convert list of strings to chars for parallelization
-                    try:
-                        unicode_type = len(max(var_dict['data'].flatten(), key=len))
-                        if ((var_dict['data'].dtype == np.dtype('<U' + str(unicode_type))) or 
-                            (var_dict['data'].dtype == np.str)):
-                            # Add strlen as a dimension if needed
-                            if self.strlen is not None:
-                                if 'strlen' not in var_dims:
-                                    var_dims += ('strlen',)
-                                var_dict['data'] = stringtochar(np.array([v.encode('ascii', 'ignore') 
-                                                                for v in var_dict['data']]).astype('S' + 
-                                                                                                   str(self.strlen)))
-                                var_dtype = 'S1'
-                    except:
-                        pass
+                    if (var_dict['data'].dtype != float) and (var_dict['data'].dtype != int):
+                        try:
+                            # Get unicode
+                            unicode_type = len(max(var_dict['data'].flatten(), key=len))
+                            
+                            if ((var_dict['data'].dtype == np.dtype('<U' + str(unicode_type))) or 
+                                (var_dict['data'].dtype == np.str)):
+                                
+                                if self.strlen is not None:
+
+                                    # Add strlen as a dimension if needed
+                                    if 'strlen' not in var_dims:
+                                        var_dims += ('strlen',)
+                                    
+                                    # Split strings into chars (S1)
+                                    var_dict['data'] = stringtochar(np.array([v.encode('ascii', 'ignore') 
+                                                                    for v in var_dict['data']]).astype('S' + 
+                                                                                                    str(self.strlen)))
+                                    var_dtype = 'S1'
+
+                        # TODO: Detect exception
+                        except Exception as e:
+                            print('String data of {0} could not be converted into chars while writing'.format(var_name))
+                            print('Error:' + e)
+                            # TODO document that exception
+                            pass
                     
                     if self.info:
                         print("Rank {0:03d}: Writing {1} var ({2}/{3})".format(self.rank, var_name, i + 1,
diff --git a/tests/4.2-test_sum.py b/tests/4.2-test_sum.py
new file mode 100644
index 0000000000000000000000000000000000000000..2fb9129e466ab838670838e040ba8eb19bea3a23
--- /dev/null
+++ b/tests/4.2-test_sum.py
@@ -0,0 +1,76 @@
+#!/usr/bin/env python
+
+import sys
+from mpi4py import MPI
+import pandas as pd
+import timeit
+import numpy as np
+from nes import *
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+size = comm.Get_size()
+
+parallel_method = 'Y'
+
+result_path = "Times_test_4.2_sum_{0}_{1:03d}.csv".format(parallel_method, size)
+result = pd.DataFrame(index=['read', 'calculate', 'write'],
+                      columns=['4.2.1.Sum'])
+
+# ======================================================================================================================
+# ===================================     CENTROID FROM NEW FILE     ===================================================
+# ======================================================================================================================
+
+test_name = '4.2.1.Sum'
+
+if rank == 0:
+    print(test_name)
+
+# CREATE GRID
+st_time = timeit.default_timer()
+projection = 'regular'
+lat_orig = 41.1
+lon_orig = 1.8
+inc_lat = 0.2
+inc_lon = 0.2
+n_lat = 100
+n_lon = 100
+nessy = create_nes(projection=projection, lat_orig=lat_orig, lon_orig=lon_orig, inc_lat=inc_lat, inc_lon=inc_lon,
+                   n_lat=n_lat, n_lon=n_lon)
+
+# ADD VARIABLES
+nessy.variables = {'var_aux': {'data': np.ones((len(nessy.time), len(nessy.lev['data']),
+                                                nessy.lat['data'].shape[0], nessy.lon['data'].shape[-1]))}}
+
+# CREATE GRID WITH COPY
+nessy_2 = nessy.copy(copy_vars=True)
+
+# ADD VARIABLES
+for var_name in nessy_2.variables.keys():
+    nessy_2.variables[var_name]['data'] *= 2
+
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+# SUM
+st_time = timeit.default_timer()
+nessy_3 = nessy + nessy_2
+print('Sum result', nessy_3.variables['var_aux']['data'])
+
+comm.Barrier()
+result.loc['calculate', test_name] = timeit.default_timer() - st_time
+
+# WRITE
+st_time = timeit.default_timer()
+nessy_3.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size), info=True)
+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/run_scalability_tests_nord3v2.sh b/tests/run_scalability_tests_nord3v2.sh
index 4e6a8fe63f92342b78e262d19dff9e70cf68cdd1..c75e60129a7c10840d594d8346052701002bb966 100644
--- a/tests/run_scalability_tests_nord3v2.sh
+++ b/tests/run_scalability_tests_nord3v2.sh
@@ -8,7 +8,7 @@ module load Python/3.7.4-GCCcore-8.3.0
 module load NES/1.0.0-nord3-v2-foss-2019b-Python-3.7.4
 
 
-for EXE in "1.1-test_read_write_projection.py" "1.2-test_create_projection.py" "1.3-test_selecting.py" "2.1-test_spatial_join.py" "2.2-test_create_shapefile.py" "2.3-test_bounds.py" "2.4-test_cell_area.py" "3.1-test_vertical_interp.py" "3.2-test_horiz_interp_bilinear.py" "3.3-test_horiz_interp_conservative.py" "4.1-test_daily_stats.py"
+for EXE in "1.1-test_read_write_projection.py" "1.2-test_create_projection.py" "1.3-test_selecting.py" "2.1-test_spatial_join.py" "2.2-test_create_shapefile.py" "2.3-test_bounds.py" "2.4-test_cell_area.py" "3.1-test_vertical_interp.py" "3.2-test_horiz_interp_bilinear.py" "3.3-test_horiz_interp_conservative.py" "4.1-test_daily_stats.py" "4.2-test_sum.py"
     do
         for nprocs in 1 2 4 8 16
             do
diff --git a/tests/test_bash_mn4.cmd b/tests/test_bash_mn4.cmd
index f8663a1d6f5d79192ab77f1418acb5e596cbe98c..491076a8110a91f178352d4f1d1d97b5d5ecad1a 100644
--- a/tests/test_bash_mn4.cmd
+++ b/tests/test_bash_mn4.cmd
@@ -18,7 +18,7 @@ module use /gpfs/projects/bsc32/software/suselinux/11/modules/all
 module load NES/1.1.0-mn4-foss-2019b-Python-3.7.4
 module load OpenMPI/4.0.5-GCC-8.3.0-mn4
 
-cd /gpfs/projects/bsc32/models/NES_master/tests
+cd /gpfs/projects/bsc32/models/NES_master/tests || exit
 
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 1.1-test_read_write_projection.py
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 1.2-test_create_projection.py
@@ -34,3 +34,4 @@ mpirun --mca mpi_warn_on_fork 0 -np 4 python 3.2-test_horiz_interp_bilinear.py
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 3.3-test_horiz_interp_conservative.py
 
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 4.1-test_daily_stats.py
+mpirun --mca mpi_warn_on_fork 0 -np 4 python 4.2-test_sum.py
diff --git a/tests/test_bash_nord3v2.cmd b/tests/test_bash_nord3v2.cmd
index 15f22c9d30942e23ebb094f55a508c8944f70c45..45a2e3db91d271b3d01a88605ebf5a0054f7ca6b 100644
--- a/tests/test_bash_nord3v2.cmd
+++ b/tests/test_bash_nord3v2.cmd
@@ -16,7 +16,7 @@ module purge
 
 module load NES/1.1.0-nord3-v2-foss-2019b-Python-3.7.4
 
-cd /gpfs/projects/bsc32/models/NES_master/tests
+cd /gpfs/projects/bsc32/models/NES_master/tests || exit
 
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 1.1-test_read_write_projection.py
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 1.2-test_create_projection.py
@@ -32,3 +32,4 @@ mpirun --mca mpi_warn_on_fork 0 -np 4 python 3.2-test_horiz_interp_bilinear.py
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 3.3-test_horiz_interp_conservative.py
 
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 4.1-test_daily_stats.py
+mpirun --mca mpi_warn_on_fork 0 -np 4 python 4.2-test_sum.py
diff --git a/tutorials/3.Statistics/3.2.Sum.ipynb b/tutorials/3.Statistics/3.2.Sum.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..53bf070bf56d4abd307d5d49c34d5b4f9c1ac754
--- /dev/null
+++ b/tutorials/3.Statistics/3.2.Sum.ipynb
@@ -0,0 +1,154 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Sum"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from nes import *\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Define first grid"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "projection = 'regular'\n",
+    "lat_orig = 41.1\n",
+    "lon_orig = 1.8\n",
+    "inc_lat = 0.2\n",
+    "inc_lon = 0.2\n",
+    "n_lat = 100\n",
+    "n_lon = 100\n",
+    "nessy = create_nes(projection=projection, lat_orig=lat_orig, lon_orig=lon_orig, inc_lat=inc_lat, inc_lon=inc_lon,\n",
+    "                   n_lat=n_lat, n_lon=n_lon)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nessy.variables = {'var_aux': {'data': np.ones((len(nessy.time), len(nessy.lev['data']),\n",
+    "                                                nessy.lat['data'].shape[0], nessy.lon['data'].shape[-1]))}}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nessy.variables"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Define second grid"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nessy_2 = nessy.copy(copy_vars=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for var_name in nessy_2.variables.keys():\n",
+    "    nessy_2.variables[var_name]['data'] *= 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nessy_2.variables"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Sum variables"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nessy_3 = nessy + nessy_2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nessy_3.variables"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nessy_3.to_netcdf(\"sum.nc\", info=True)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}