diff --git a/nes/load_nes.py b/nes/load_nes.py
index 5bfd333fc4661ec1f9feffe82dfe4d726051b13e..8de2346d2ccce85d61048f18df922f3f5e959986 100644
--- a/nes/load_nes.py
+++ b/nes/load_nes.py
@@ -4,9 +4,12 @@ import os
 from mpi4py import MPI
 from netCDF4 import Dataset
 import warnings
-
+import numpy as np
 from .nc_projections import *
 
+DIM_VAR_NAMES = ['lat', 'latitude', 'lat_bnds', 'lon', 'longitude', 'lon_bnds', 'time', 'time_bnds', 'lev', 'level',
+                 'cell_area', 'crs', 'rotated_pole', 'x', 'y', 'rlat', 'rlon', 'Lambert_conformal', 'mercator']
+
 
 def open_netcdf(path, comm=None, xarray=False, info=False, parallel_method='Y', avoid_first_hours=0, avoid_last_hours=0,
                 first_level=0, last_level=None, balanced=False):
@@ -270,15 +273,53 @@ def concatenate_netcdfs(nessy_list, comm=None, info=False, parallel_method='Y',
         nessy_first = nessy_list[0]
     for i, aux_nessy in enumerate(nessy_list[1:]):
         if isinstance(aux_nessy, str):
-            aux_nessy = open_netcdf(aux_nessy,
-                                    comm=comm,
-                                    parallel_method=parallel_method,
-                                    avoid_first_hours=avoid_first_hours,
-                                    avoid_last_hours=avoid_last_hours,
-                                    first_level=first_level,
-                                    last_level=last_level,
-                                    balanced=balanced
-                                    )
-        nessy_first.concatenate(aux_nessy)
+            nc_add = Dataset(filename=aux_nessy, mode='r')
+            for var_name, var_info in nc_add.variables.items():
+                if var_name not in DIM_VAR_NAMES:
+                    nessy_first.variables[var_name] = {}
+                    var_dims = var_info.dimensions
+                    # Read data in 4 dimensions
+                    if len(var_dims) < 2:
+                        data = var_info[:]
+                    elif len(var_dims) == 2:
+                        data = var_info[nessy_first.read_axis_limits['y_min']:nessy_first.read_axis_limits['y_max'],
+                                        nessy_first.read_axis_limits['x_min']:nessy_first.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 = var_info[nessy_first.read_axis_limits['y_min']:nessy_first.read_axis_limits['y_max'],
+                                            nessy_first.read_axis_limits['x_min']:nessy_first.read_axis_limits['x_max'],
+                                            :]
+                            data_aux = np.empty(shape=(data.shape[0], data.shape[1]), dtype=np.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].tostring().decode('ascii').replace('\x00', ''))
+                            data = data_aux.reshape((1, 1, data_aux.shape[-2], data_aux.shape[-1]))
+                        else:
+                            data = var_info[nessy_first.read_axis_limits['t_min']:nessy_first.read_axis_limits['t_max'],
+                                            nessy_first.read_axis_limits['y_min']:nessy_first.read_axis_limits['y_max'],
+                                            nessy_first.read_axis_limits['x_min']:nessy_first.read_axis_limits['x_max']]
+                            data = data.reshape(data.shape[-3], 1, data.shape[-2], data.shape[-1])
+                    elif len(var_dims) == 4:
+                        data = var_info[nessy_first.read_axis_limits['t_min']:nessy_first.read_axis_limits['t_max'],
+                                        nessy_first.read_axis_limits['z_min']:nessy_first.read_axis_limits['z_max'],
+                                        nessy_first.read_axis_limits['y_min']:nessy_first.read_axis_limits['y_max'],
+                                        nessy_first.read_axis_limits['x_min']:nessy_first.read_axis_limits['x_max']]
+                    else:
+                        raise TypeError("{} data shape is nto accepted".format(var_dims))
+
+                    nessy_first.variables[var_name]['data'] = data
+                    # Avoid some attributes
+                    for attrname in var_info.ncattrs():
+                        if attrname not in ['missing_value', '_FillValue']:
+                            value = getattr(var_info, attrname)
+                            if value in ['unitless', '-']:
+                                value = ''
+                            nessy_first.variables[var_name][attrname] = value
+            nc_add.close()
+
+        else:
+            nessy_first.concatenate(aux_nessy)
         
     return nessy_first
diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py
index 2d056372d0a2a88377f77bd0be6b94b870d3ef24..9585c790cfd3f5b452697e36000ab9baa6375c42 100644
--- a/nes/nc_projections/default_nes.py
+++ b/nes/nc_projections/default_nes.py
@@ -2238,66 +2238,68 @@ class Nes(object):
 
         for i, (var_name, var_dict) in enumerate(self.variables.items()):
             if var_dict['data'] is not None:
-                
-                # Get dimensions
-                if len(var_dict['data'].shape) == 4:
+                if isinstance(var_dict['data'], int) and var_dict['data'] == 0:
                     var_dims = ('time', 'lev',) + self._var_dim
+                    var_dtype = np.float32
                 else:
-                    var_dims = self._var_dim
-
-                # Get data type
-                if 'dtype' in var_dict.keys():
-                    var_dtype = var_dict['dtype']
-                    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)
-                        warnings.warn(msg)
+                    # Get dimensions
+                    if len(var_dict['data'].shape) == 4:
+                        var_dims = ('time', 'lev',) + self._var_dim
+                    else:
+                        var_dims = self._var_dim
+
+                    # Get data type
+                    if 'dtype' in var_dict.keys():
+                        var_dtype = var_dict['dtype']
+                        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)
+                            warnings.warn(msg)
+                            try:
+                                var_dict['data'] = var_dict['data'].astype(var_dtype)
+                            except Exception as e:  # TODO: Detect exception
+                                print(e)
+                                raise TypeError("It was not possible to cast the data to the input dtype.")
+                    else:
+                        var_dtype = var_dict['data'].dtype
+                    # Transform objects into strings
+                    if var_dtype == np.dtype(object):
+                        var_dict['data'] = var_dict['data'].astype(str)
+                        var_dtype = var_dict['data'].dtype
+
+                    # Convert list of strings to chars for parallelization
+                    if not np.issubdtype(var_dict['data'].dtype, np.number):
                         try:
-                            var_dict['data'] = var_dict['data'].astype(var_dtype)
-                        except Exception as e:  # TODO: Detect exception
+                            # 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 == 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(e)
-                            raise TypeError("It was not possible to cast the data to the input dtype.")
-                else:
-                    var_dtype = var_dict['data'].dtype
-
-                # Transform objects into strings
-                if var_dtype == np.dtype(object):
-                    var_dict['data'] = var_dict['data'].astype(str)
-                    var_dtype = var_dict['data'].dtype
-
-                # Convert list of strings to chars for parallelization
-                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 == 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
+                            # TODO document that exception
+                            pass
 
                 if self.info:
                     print("Rank {0:03d}: Writing {1} var ({2}/{3})".format(
@@ -2330,7 +2332,12 @@ class Nes(object):
                         if att_name == 'data':
                             if self.info:
                                 print("Rank {0:03d}: Filling {1})".format(self.rank, var_name))
-                            if len(att_value.shape) == 4:
+                            if isinstance(att_value, int) and att_value == 0:
+                                var[self.write_axis_limits['t_min']:self.write_axis_limits['t_max'],
+                                    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']] = 0
+                            elif len(att_value.shape) == 4:
                                 try:
                                     var[self.write_axis_limits['t_min']:self.write_axis_limits['t_max'],
                                         self.write_axis_limits['z_min']:self.write_axis_limits['z_max'],
@@ -3067,8 +3074,8 @@ class Nes(object):
         for var_name in data_list.keys():
             if self.info and self.master:
                 print("Gathering {0}".format(var_name))
-            shp_len = len(data_list[var_name]['data'].shape)
             try:
+                shp_len = len(data_list[var_name]['data'].shape)
                 # 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()
@@ -3127,6 +3134,9 @@ class Nes(object):
                         data_list[var_name]['data'] = np.stack(recvbuf)
                     else:
                         data_list[var_name]['data'] = np.concatenate(recvbuf, axis=axis)
+
+            except AttributeError:
+                data_list[var_name]['data'] = 0
             except Exception as e:
                 print("**ERROR** an error has occurred while gathering the '{0}' variable.\n".format(var_name))
                 sys.stderr.write("**ERROR** an error has occurred while gathering the '{0}' variable.\n".format(
diff --git a/nes/nc_projections/points_nes.py b/nes/nc_projections/points_nes.py
index 318302e23519ac3ce5ae0fb9e58353fc8a0b0de2..f542bacf94b487d990fc198b7fabb8286f9f1ac3 100644
--- a/nes/nc_projections/points_nes.py
+++ b/nes/nc_projections/points_nes.py
@@ -368,7 +368,7 @@ class PointsNes(Nes):
                             var_dims = ('time',) + self._var_dim
 
                     # Convert list of strings to chars for parallelization
-                    if (var_dict['data'].dtype != float) and (var_dict['data'].dtype != int):
+                    if not np.issubdtype(var_dict['data'].dtype, np.number):
                         try:
                             # Get unicode 
                             unicode_type = len(max(var_dict['data'].flatten(), key=len))
diff --git a/nes/nc_projections/points_nes_ghost.py b/nes/nc_projections/points_nes_ghost.py
index 695ac4b658dc36b38f9434f980b470eafc2e2ff6..936cbca5346fbf0f5d79c2d21a893fea710b1901 100644
--- a/nes/nc_projections/points_nes_ghost.py
+++ b/nes/nc_projections/points_nes_ghost.py
@@ -350,7 +350,7 @@ class PointsNesGHOST(PointsNes):
                         var_dims = self._var_dim + ('time',)
 
                     # Convert list of strings to chars for parallelization
-                    if (var_dict['data'].dtype != float) and (var_dict['data'].dtype != int):
+                    if not np.issubdtype(var_dict['data'].dtype, np.number):
                         try:
                             # Get unicode
                             unicode_type = len(max(var_dict['data'].flatten(), key=len))
diff --git a/nes/nc_projections/points_nes_providentia.py b/nes/nc_projections/points_nes_providentia.py
index 3e45308b688876d935682f1cfff3ac8da66e02b0..937d7e676211d5837e99254aebcfebf57209f12b 100644
--- a/nes/nc_projections/points_nes_providentia.py
+++ b/nes/nc_projections/points_nes_providentia.py
@@ -386,7 +386,7 @@ class PointsNesProvidentia(PointsNes):
                         var_dims = self._var_dim + ('time',)
 
                     # Convert list of strings to chars for parallelization
-                    if (var_dict['data'].dtype != float) and (var_dict['data'].dtype != int):
+                    if not np.issubdtype(var_dict['data'].dtype, np.number):
                         try:
                             # Get unicode
                             unicode_type = len(max(var_dict['data'].flatten(), key=len))