diff --git a/nes/methods/horizontal_interpolation.py b/nes/methods/horizontal_interpolation.py
index 8f480345426e6ef8e3bc64a32a763125356264c5..5fc29d87aa68f277eb18744349feada904e72291 100644
--- a/nes/methods/horizontal_interpolation.py
+++ b/nes/methods/horizontal_interpolation.py
@@ -3,6 +3,7 @@
 import sys
 import warnings
 import numpy as np
+import pandas as pd
 import os
 import nes
 from mpi4py import MPI
@@ -10,10 +11,16 @@ from scipy import spatial
 from filelock import FileLock
 from datetime import datetime
 from warnings import warn
+import copy
+import pyproj
+
+# CONSTANTS
+NEAREST_OPTS = ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']
+CONSERVATIVE_OPTS = ['Conservative', 'Area_Conservative', 'cons', 'conservative', 'area']
 
 
 def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='NearestNeighbour', n_neighbours=4,
-                           info=False, to_providentia=False):
+                           info=False, to_providentia=False, only_create_wm=False, wm=None):
     """
     Horizontal methods from one grid to another one.
 
@@ -26,23 +33,44 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares
     weight_matrix_path : str, None
         Path to the weight matrix to read/create.
     kind : str
-        Kind of horizontal methods. Accepted values: ['NearestNeighbour'].
+        Kind of horizontal methods. Accepted values: ['NearestNeighbour', 'Conservative'].
     n_neighbours : int
         Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4.
     info : bool
         Indicates if you want to print extra info during the methods process.
     to_providentia : bool
         Indicates if we want the interpolated grid in Providentia format.
+    only_create_wm : bool
+        Indicates if you want to only create the Weight Matrix.
+    wm : Nes
+        Weight matrix Nes File
     """
-
+    if info and self.master:
+        print("Creating Weight Matrix")
     # Obtain weight matrix
     if self.parallel_method == 'T':
-        weights, idx = get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours)
+        weights, idx = get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours,
+                                              only_create_wm, wm)
     elif self.parallel_method in ['Y', 'X']:
-        weights, idx = get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours)
+        weights, idx = get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours,
+                                               only_create_wm, wm)
     else:
         raise NotImplemented("Parallel method {0} is not implemented yet for horizontal interpolations. Use 'T'".format(
             self.parallel_method))
+    if info and self.master:
+        print("Weight Matrix done!")
+    if only_create_wm:
+        # weights for only_create is the WM NES object
+        return weights
+
+    # idx[idx < 0] = np.nan
+    idx = np.ma.masked_array(idx, mask=idx == -999)
+    # idx = np.array(idx, dtype=float)
+    # idx[idx < 0] = np.nan
+    # weights[weights < 0] = np.nan
+    weights = np.ma.masked_array(weights, mask=weights == -999)
+    # weights = np.array(weights, dtype=float)
+    # weights[weights < 0] = np.nan
 
     # Copy NES
     final_dst = dst_grid.copy()
@@ -61,6 +89,8 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares
     final_dst.hours_start = self.hours_start
     final_dst.hours_end = self.hours_end
 
+    if info and self.master:
+        print("Applying weights")
     # Apply weights
     for var_name, var_info in self.variables.items():
         if info and self.master:
@@ -81,8 +111,8 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares
         for time in range(dst_shape[0]):
             for lev in range(dst_shape[1]):
                 src_aux = get_src_data(self.comm, var_info['data'][time, lev], idx, self.parallel_method)
-                # src_aux = np.take(src_data[time, lev], idx)
                 final_dst.variables[var_name]['data'][time, lev] = np.sum(weights * src_aux, axis=1)
+
         if isinstance(dst_grid, nes.PointsNes):
             # Removing level axis
             if src_shape[1] != 1:
@@ -94,6 +124,9 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares
                 
     final_dst.global_attrs = self.global_attrs
 
+    if info and self.master:
+        print("Formatting")
+
     if to_providentia:
         # self = experiment to interpolate (regular, rotated, etc.)
         # final_dst = interpolated experiment (points)
@@ -154,13 +187,14 @@ def get_src_data(comm, var_data, idx, parallel_method):
 
         var_data = comm.bcast(var_data)
 
-    var_data = np.take(var_data, idx)
+    var_data = np.pad(var_data, [1, 1], 'constant', constant_values=np.nan).take(idx + 1, mode='clip')
+    #var_data = np.take(var_data, idx)
 
     return var_data
 
 
 # noinspection DuplicatedCode
-def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours):
+def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours, only_create, wm):
     """
     To obtain the weights and source data index through the T axis.
 
@@ -173,51 +207,75 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour
     weight_matrix_path : str, None
         Path to the weight matrix to read/create.
     kind : str
-        Kind of horizontal methods. Accepted values: ['NearestNeighbour'].
+        Kind of horizontal methods. Accepted values: ['NearestNeighbour', 'Conservative'].
     n_neighbours : int
         Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4.
+    only_create : bool
+        Indicates if you want to only create the Weight Matrix.
+    wm : Nes
+        Weight matrix Nes File
 
     Returns
     -------
     tuple
         Weights and source data index.
     """
+    if wm is not None:
+        weight_matrix = wm
 
-    if weight_matrix_path is not None:
-        with FileLock(weight_matrix_path.replace('.nc', '.lock')):
-            if self.master:
-                if os.path.isfile(weight_matrix_path):
+    elif weight_matrix_path is not None:
+        with FileLock(weight_matrix_path + "{0:03d}.lock".format(self.rank)):
+            if os.path.isfile(weight_matrix_path):
+                if self.master:
                     weight_matrix = read_weight_matrix(weight_matrix_path, comm=MPI.COMM_SELF)
-                    if len(weight_matrix.lev['data']) != n_neighbours:
-                        warn("The selected weight matrix does not have the same number of nearest neighbours." +
-                             "Re-calculating again but not saving it.")
-                        if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']:
-                            weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours)
-                        else:
-                            raise NotImplementedError(kind)
                 else:
-                    if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']:
-                        weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours)
+                    weight_matrix = True
+                if kind in NEAREST_OPTS:
+                    if self.master:
+                        if len(weight_matrix.lev['data']) != n_neighbours:
+                            warn("The selected weight matrix does not have the same number of nearest neighbours." +
+                                 "Re-calculating again but not saving it.")
+                            weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours)
                     else:
-                        raise NotImplementedError(kind)
-                    if weight_matrix_path is not None:
-                        weight_matrix.to_netcdf(weight_matrix_path)
+                        weight_matrix = True
+
             else:
-                weight_matrix = None
+                if self.master:
+                    if kind in NEAREST_OPTS:
+                        weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours,
+                                                                wm_path=weight_matrix_path)
+                    elif kind in CONSERVATIVE_OPTS:
+                        weight_matrix = create_area_conservative_weight_matrix(self, dst_grid,
+                                                                               wm_path=weight_matrix_path)
+                    else:
+                        raise NotImplementedError(kind)
+                else:
+                    weight_matrix = True
+
+        if os.path.exists(weight_matrix_path + "{0:03d}.lock".format(self.rank)):
+            os.remove(weight_matrix_path + "{0:03d}.lock".format(self.rank))
     else:
-        if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']:
-            if self.master:
+        if self.master:
+            if kind in NEAREST_OPTS:
                 weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours)
+            elif kind in CONSERVATIVE_OPTS:
+                weight_matrix = create_area_conservative_weight_matrix(self, dst_grid)
             else:
-                weight_matrix = None
+                raise NotImplementedError(kind)
         else:
-            raise NotImplementedError(kind)
+            weight_matrix = True
+
+    if only_create:
+        return weight_matrix, None
 
-    # Normalize to 1
     if self.master:
-        weights = np.array(np.array(weight_matrix.variables['inverse_dists']['data'], dtype=np.float64) /
-                           np.array(weight_matrix.variables['inverse_dists']['data'], dtype=np.float64).sum(axis=1),
-                           dtype=np.float64)
+        if kind in NEAREST_OPTS:
+            # Normalize to 1
+            weights = np.array(np.array(weight_matrix.variables['weight']['data'], dtype=np.float64) /
+                               np.array(weight_matrix.variables['weight']['data'], dtype=np.float64).sum(axis=1),
+                               dtype=np.float64)
+        else:
+            weights = np.array(weight_matrix.variables['weight']['data'], dtype=np.float64)
         idx = np.array(weight_matrix.variables['idx']['data'][0], dtype=int)
     else:
         weights = None
@@ -230,7 +288,7 @@ def get_weights_idx_t_axis(self, dst_grid, weight_matrix_path, kind, n_neighbour
 
 
 # noinspection DuplicatedCode
-def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours):
+def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbours, only_create, wm):
     """
     To obtain the weights and source data index through the X or Y axis.
 
@@ -243,63 +301,83 @@ def get_weights_idx_xy_axis(self, dst_grid, weight_matrix_path, kind, n_neighbou
     weight_matrix_path : str, None
         Path to the weight matrix to read/create.
     kind : str
-        Kind of horizontal methods. Accepted values: ['NearestNeighbour']
+        Kind of horizontal methods. Accepted values: ['NearestNeighbour', 'Conservative'].
     n_neighbours : int
         Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4.
+    only_create : bool
+        Indicates if you want to only create the Weight Matrix.
+    wm : Nes
+        Weight matrix Nes File
 
     Returns
     -------
     tuple
         Weights and source data index.
     """
-
     if isinstance(dst_grid, nes.PointsNes) and weight_matrix_path is not None:
         if self.master:
             warn("To point weight matrix cannot be saved.")
         weight_matrix_path = None
 
-    if weight_matrix_path is not None:
-        with FileLock(weight_matrix_path.replace('.nc', '.lock')):
-            if self.master:
-                if os.path.isfile(weight_matrix_path):
+    if wm is not None:
+       weight_matrix = wm
+
+    elif weight_matrix_path is not None:
+        with FileLock(weight_matrix_path + "{0:03d}.lock".format(self.rank)):
+            if os.path.isfile(weight_matrix_path):
+                if self.master:
                     weight_matrix = read_weight_matrix(weight_matrix_path, comm=MPI.COMM_SELF)
+                else:
+                    weight_matrix = True
+                if kind in NEAREST_OPTS:
                     if len(weight_matrix.lev['data']) != n_neighbours:
                         warn("The selected weight matrix does not have the same number of nearest neighbours." +
                              "Re-calculating again but not saving it.")
-                        if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']:
+                        if self.master:
                             weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours)
                         else:
-                            raise NotImplementedError(kind)
-                else:
-                    if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']:
-                        weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours)
-                    else:
-                        raise NotImplementedError(kind)
-                    if weight_matrix_path is not None:
-                        weight_matrix.to_netcdf(weight_matrix_path)
+                            weight_matrix = True
             else:
-                weight_matrix = None
+                if kind in NEAREST_OPTS:
+                    if self.master:
+                        weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours,
+                                                                wm_path=weight_matrix_path)
+                    else:
+                        weight_matrix = True
+                elif kind in CONSERVATIVE_OPTS:
+                    weight_matrix = create_area_conservative_weight_matrix(self, dst_grid, wm_path=weight_matrix_path)
+                else:
+                    raise NotImplementedError(kind)
+
+        if os.path.exists(weight_matrix_path + "{0:03d}.lock".format(self.rank)):
+            os.remove(weight_matrix_path + "{0:03d}.lock".format(self.rank))
     else:
-        if kind in ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN']:
-            if self.master:
-                weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours)
-            else:
-                weight_matrix = None
+        if kind in NEAREST_OPTS:
+            weight_matrix = create_nn_weight_matrix(self, dst_grid, n_neighbours=n_neighbours)
+        elif kind in CONSERVATIVE_OPTS:
+            weight_matrix = create_area_conservative_weight_matrix(self, dst_grid)
         else:
             raise NotImplementedError(kind)
 
+    if only_create:
+        return weight_matrix, None
+
     # Normalize to 1
     if self.master:
-        weights = np.array(np.array(weight_matrix.variables['inverse_dists']['data'], dtype=np.float64) /
-                           np.array(weight_matrix.variables['inverse_dists']['data'], dtype=np.float64).sum(axis=1),
-                           dtype=np.float64)
-        idx = np.array(weight_matrix.variables['idx']['data'][0], dtype=int)
+        if kind in NEAREST_OPTS:
+            weights = np.array(np.array(weight_matrix.variables['weight']['data'], dtype=np.float64) /
+                               np.array(weight_matrix.variables['weight']['data'], dtype=np.float64).sum(axis=1),
+                               dtype=np.float64)
+        else:
+            weights = np.array(weight_matrix.variables['weight']['data'], dtype=np.float64)
+        idx = np.array(weight_matrix.variables['idx']['data'][0], dtype=np.int64)
     else:
         weights = None
         idx = None
 
     weights = self.comm.bcast(weights, root=0)
     idx = self.comm.bcast(idx, root=0)
+
     # if isinstance(dst_grid, nes.PointsNes):
     # print("weights 1 ->", weights.shape)
     # print("idx 1 ->", idx.shape)
@@ -334,14 +412,17 @@ def read_weight_matrix(weight_matrix_path, comm=None, parallel_method='T'):
     nes.Nes
         Weight matrix.
     """
-
     weight_matrix = nes.open_netcdf(path=weight_matrix_path, comm=comm, parallel_method=parallel_method, balanced=True)
     weight_matrix.load()
 
+    # In previous versions of NES weight was called inverse_dists
+    if 'inverse_dists' in weight_matrix.variables.keys():
+        weight_matrix.variables['weight'] = weight_matrix.variables['inverse_dists']
+
     return weight_matrix
 
 
-def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False):
+def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, wm_path=None, info=False):
     """
     To create the weight matrix with the nearest neighbours method.
 
@@ -390,8 +471,12 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False):
     # model geographic longitude/latitude coordinates are first converted
     # to cartesian ECEF (Earth Centred, Earth Fixed) coordinates, before
     # calculating distances.
-    src_mod_xy = self.lon_lat_to_cartesian_ecef(src_lon.flatten(), src_lat.flatten())
-    dst_mod_xy = self.lon_lat_to_cartesian_ecef(dst_lon.flatten(), dst_lat.flatten())
+
+    # src_mod_xy = lon_lat_to_cartesian(src_lon.flatten(), src_lat.flatten())
+    # dst_mod_xy = lon_lat_to_cartesian(dst_lon.flatten(), dst_lat.flatten())
+
+    src_mod_xy = lon_lat_to_cartesian_ecef(src_lon.flatten(), src_lat.flatten())
+    dst_mod_xy = lon_lat_to_cartesian_ecef(dst_lon.flatten(), dst_lat.flatten())
 
     # generate KDtree using model 1 coordinates (i.e. the model grid you are
     # interpolating from)
@@ -407,6 +492,8 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False):
     weight_matrix = dst_grid.copy()
     weight_matrix.time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)]
     weight_matrix._time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)]
+    weight_matrix._time_bnds = None
+    weight_matrix.time_bnds = None
     weight_matrix.last_level = None
     weight_matrix.first_level = 0
     weight_matrix.hours_start = 0
@@ -418,10 +505,216 @@ def create_nn_weight_matrix(self, dst_grid, n_neighbours=4, info=False):
     inverse_dists = np.reciprocal(dists)
 
     inverse_dists_transf = inverse_dists.T.reshape((1, n_neighbours, dst_lon.shape[0], dst_lon.shape[1]))
-    weight_matrix.variables['inverse_dists'] = {'data': inverse_dists_transf, 'units': 'm'}
+    weight_matrix.variables['weight'] = {'data': inverse_dists_transf, 'units': 'm'}
     idx_transf = idx.T.reshape((1, n_neighbours, dst_lon.shape[0], dst_lon.shape[1]))
     weight_matrix.variables['idx'] = {'data': idx_transf, 'units': ''}
     weight_matrix.lev = {'data': np.arange(inverse_dists_transf.shape[1]), 'units': ''}
     weight_matrix._lev = {'data': np.arange(inverse_dists_transf.shape[1]), 'units': ''}
+    if wm_path is not None:
+        weight_matrix.to_netcdf(wm_path)
+    
+    return weight_matrix
+
+
+def create_area_conservative_weight_matrix(self, dst_nes, wm_path=None, info=False):
+    """
+    To create the weight matrix with the area conservative method.
+
+    Parameters
+    ----------
+    self : nes.Nes
+        Source projection Nes Object.
+    dst_nes : nes.Nes
+        Final projection Nes object.
+
+    info: bool
+        Indicates if you want to print extra info during the methods process.
+
+    Returns
+    -------
+    nes.Nes
+        Weight matrix.
+    """
+    if info and self.master:
+        print("\tCreating area conservative Weight Matrix")
+        sys.stdout.flush()
+    # Get a portion of the destiny grid
+    if dst_nes.shapefile is None:
+        dst_nes.create_shapefile()
+    dst_grid = copy.deepcopy(dst_nes.shapefile)
+
+    # Get the complete source grid
+    if self.shapefile is None:
+        self.create_shapefile()
+    src_grid = copy.deepcopy(self.shapefile)
+
+    if self.parallel_method == 'T':
+        # All process has the same shapefile
+        pass
+    else:
+        src_grid = self.comm.gather(src_grid, root=0)
+        if self.master:
+            src_grid = pd.concat(src_grid)
+        src_grid = self.comm.bcast(src_grid)
+
+    my_crs = pyproj.CRS.from_proj4("+proj=latlon")
+    # Normalizing projections
+    dst_grid.to_crs(crs=my_crs, inplace=True)
+    dst_grid['FID_dst'] = dst_grid.index
+    dst_grid = dst_grid.reset_index()
+
+    # src_grid.to_crs(crs=pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84').crs, inplace=True)
+    src_grid.to_crs(crs=my_crs, inplace=True)
+    src_grid['FID_src'] = src_grid.index
+
+    src_grid = src_grid.reset_index()
+
+    if info and self.master:
+        print("\t\tGrids created and ready to interpolate")
+        sys.stdout.flush()
+
+    # Get intersected areas
+    inp, res = dst_grid.sindex.query_bulk(src_grid.geometry, predicate='intersects')
+
+    # Calculate intersected areas and fractions
+    intersection = pd.DataFrame(columns=["FID_src", "FID_dst"])
+    intersection['INP'] = np.array(inp, dtype=np.uint32)
+    intersection['RES'] = np.array(res, dtype=np.uint32)
+    intersection['FID_src'] = np.array(src_grid.loc[inp, 'FID_src'], dtype=np.uint32)
+    intersection['FID_dst'] = np.array(dst_grid.loc[res, 'FID_dst'], dtype=np.uint32)
+
+    intersection['geometry_src'] = src_grid.loc[inp, 'geometry'].values
+    intersection['geometry_dst'] = dst_grid.loc[res, 'geometry'].values
+
+    if True:
+        # No Warnings Zone
+        warnings.filterwarnings('ignore')
+        intersection['intersect_area'] = intersection.apply(
+            lambda x: x['geometry_src'].intersection(x['geometry_dst']).buffer(0).area, axis=1)
+        intersection.drop(intersection[intersection['intersect_area'] <= 0].index, inplace=True)
+
+        intersection["src_area"] = src_grid.loc[intersection['FID_src'], 'geometry'].area.values
+
+        warnings.filterwarnings('default')
+
+    intersection['weight'] = intersection['intersect_area'] / intersection["src_area"]
+
+    # Format & Clean
+    intersection.drop(columns=["geometry_src", "geometry_dst", "src_area", "intersect_area"], inplace=True)
+
+    if info and self.master:
+        print("\t\tWeights calculated. Formatting weight matrix.")
+        sys.stdout.flush()
+
+    # Initialising weight matrix
+    if self.parallel_method != 'T':
+        intersection = self.comm.gather(intersection, root=0)
+    if self.master:
+        if self.parallel_method != 'T':
+            intersection = pd.concat(intersection)
+        intersection = intersection.set_index(['FID_dst', intersection.groupby('FID_dst').cumcount()]).rename_axis(
+            ('FID', 'level')).sort_index()
+        intersection.rename(columns={"FID_src": "idx"}, inplace=True)
+        weight_matrix = dst_nes.copy()
+        weight_matrix.time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)]
+        weight_matrix._time = [datetime(year=2000, month=1, day=1, hour=0, second=0, microsecond=0)]
+        weight_matrix._time_bnds = None
+        weight_matrix.time_bnds = None
+        weight_matrix.last_level = None
+        weight_matrix.first_level = 0
+        weight_matrix.hours_start = 0
+        weight_matrix.hours_end = 0
+    
+        weight_matrix.set_communicator(MPI.COMM_SELF)
+    
+        weight_matrix.set_levels({'data': np.arange(intersection.index.get_level_values('level').max() + 1),
+                                  'dimensions': ('lev',),
+                                  'units': '',
+                                  'positive': 'up'})
+    
+        # Creating Weight matrix empty variables
+        if len(weight_matrix._lat['data'].shape) == 1:
+            shape = (1, len(weight_matrix.lev['data']),
+                     weight_matrix._lat['data'].shape[0], weight_matrix._lon['data'].shape[0],)
+            shape_flat = (1, len(weight_matrix.lev['data']),
+                          weight_matrix._lat['data'].shape[0] * weight_matrix._lon['data'].shape[0],)
+        else:
+            shape = (1, len(weight_matrix.lev['data']),
+                     weight_matrix._lat['data'].shape[0], weight_matrix._lat['data'].shape[1],)
+            shape_flat = (1, len(weight_matrix.lev['data']),
+                          weight_matrix._lat['data'].shape[0] * weight_matrix._lat['data'].shape[1],)
+    
+        weight_matrix.variables['weight'] = {'data': np.empty(shape_flat), 'units': '-'}
+        weight_matrix.variables['weight']['data'][:] = -999
+        weight_matrix.variables['idx'] = {'data': np.empty(shape_flat), 'units': '-'}
+        weight_matrix.variables['idx']['data'][:] = -999
     
+        # Filling Weight matrix variables
+        for aux_lev in weight_matrix.lev['data']:
+            aux_data = intersection.xs(level='level', key=aux_lev)
+            weight_matrix.variables['weight']['data'][0, aux_lev, aux_data.index] = aux_data.loc[:, 'weight'].values
+            weight_matrix.variables['idx']['data'][0, aux_lev, aux_data.index] = aux_data.loc[:, 'idx'].values
+        # Re-shaping
+        weight_matrix.variables['weight']['data'] = weight_matrix.variables['weight']['data'].reshape(shape)
+        weight_matrix.variables['idx']['data'] = weight_matrix.variables['idx']['data'].reshape(shape)
+        if wm_path is not None:
+            if info and self.master:
+                print("\t\tWeight matrix saved at {0}".format(wm_path))
+                sys.stdout.flush()
+            weight_matrix.to_netcdf(wm_path)
+    else:
+        weight_matrix = True
     return weight_matrix
+
+
+def lon_lat_to_cartesian(lon, lat, radius=6378137.0):
+    """
+    Calculate lon, lat coordinates of a point on a sphere.
+
+    DEPRECATED!!!!
+
+    Parameters
+    ----------
+    lon : np.array
+        Longitude values.
+    lat : np.array
+        Latitude values.
+    radius : float
+        Radius of the sphere to get the distances.
+    """
+
+    lon_r = np.radians(lon)
+    lat_r = np.radians(lat)
+
+    x = radius * np.cos(lat_r) * np.cos(lon_r)
+    y = radius * np.cos(lat_r) * np.sin(lon_r)
+    z = radius * np.sin(lat_r)
+
+    return np.column_stack([x, y, z])
+
+
+def lon_lat_to_cartesian_ecef(lon, lat):
+    """
+    # convert observational/model geographic longitude/latitude coordinates to cartesian ECEF (Earth Centred,
+    # Earth Fixed) coordinates, assuming WGS84 datum and ellipsoid, and that all heights = 0.
+    # ECEF coordiantes represent positions (in meters) as X, Y, Z coordinates, approximating the earth surface
+    # as an ellipsoid of revolution.
+    # This conversion is for the subsequent calculation of euclidean distances of the model gridcell centres
+    # from each observational station.
+    # Defining the distance between two points on the earth's surface as simply the euclidean distance
+    # between the two lat/lon pairs could lead to inaccurate results depending on the distance
+    # between two points (i.e. 1 deg. of longitude varies with latitude).
+
+    Parameters
+    ----------
+    lon : np.array
+        Longitude values.
+    lat : np.array
+        Latitude values.
+    """
+    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
+    ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
+
+    x, y, z = pyproj.transform(lla, ecef, lon, lat, np.zeros(lon.shape), radians=False)
+
+    return np.column_stack([x, y, z])
diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py
index 2686113ae3978ee1e73d548ef27f01d7caa6a7d5..1c80b9be21453d8ae0b18f4c9f22c445aaacb34c 100644
--- a/nes/nc_projections/default_nes.py
+++ b/nes/nc_projections/default_nes.py
@@ -190,7 +190,7 @@ class Nes(object):
 
             # Cell measures screening
             self.cell_measures = self.__get_cell_measures(create_nes)
-            
+
             # Set NetCDF attributes
             self.global_attrs = self.__get_global_attributes(create_nes)
 
@@ -239,7 +239,7 @@ class Nes(object):
 
             # Cell measures screening
             self.cell_measures = self._get_cell_measures_values(self._cell_measures)
-            
+
             # Set axis limits for parallel writing
             self.write_axis_limits = self.get_write_axis_limits()
 
@@ -1016,7 +1016,7 @@ class Nes(object):
                        't_min': None, 't_max': None}
         
         idx = self.get_idx_intervals()
-        
+
         if self.parallel_method == 'Y':
             y_len = idx['idx_y_max'] - idx['idx_y_min']
             if y_len < self.size:
@@ -1054,7 +1054,7 @@ class Nes(object):
             axis_limits['t_max'] = idx['idx_t_max']
             
         elif self.parallel_method == 'T':
-            t_len = idx['idx_t_min'] - idx['idx_t_max']
+            t_len = idx['idx_t_max'] - idx['idx_t_min']
             if t_len < self.size:
                 raise IndexError('More processors (size={0}) selected than T elements (size={1})'.format(
                     self.size, t_len))
@@ -1359,7 +1359,7 @@ class Nes(object):
             CF compliant time.
         """
 
-        units = time.units
+        units = self.__parse_time_unit(time.units)
         if not hasattr(time, 'calendar'):
             calendar = 'standard'
         else:
@@ -1369,7 +1369,12 @@ class Nes(object):
             units = 'days since ' + time.units.split('since')[1].lstrip()
             time = self.__get_dates_from_months(time, units, calendar)
 
-        return time, units, calendar
+        time_data = time[:]
+
+        if len(time_data) == 1 and np.isnan(time_data[0]):
+            time_data[0] = 0
+
+        return time_data, units, calendar
 
     @staticmethod
     def __parse_time_unit(t_units):
@@ -1408,8 +1413,9 @@ class Nes(object):
         else:
             if self.master:
                 nc_var = self.netcdf.variables['time']
-                nc_var, units, calendar = self.__parse_time(nc_var)
-                time = num2date(nc_var[:], self.__parse_time_unit(units), calendar=calendar)
+                time_data, units, calendar = self.__parse_time(nc_var)
+
+                time = num2date(time_data, units, calendar=calendar)
                 time = [aux.replace(second=0, microsecond=0) for aux in time]
             else:
                 time = None
@@ -1595,12 +1601,12 @@ class Nes(object):
         coordinate_len = len(values['data'].shape)
         if bounds:
             coordinate_len -= 1
-        
+
         if coordinate_axis == 'Y':
             if coordinate_len == 1:
                 values['data'] = values['data'][self.read_axis_limits['y_min']:self.read_axis_limits['y_max']]
             elif coordinate_len == 2:
-                values['data'] = values['data'][self.read_axis_limits['y_min']:self.read_axis_limits['y_max'], 
+                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(
@@ -1609,7 +1615,7 @@ class Nes(object):
             if coordinate_len == 1:
                 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'], 
+                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(
@@ -1640,25 +1646,25 @@ class Nes(object):
 
         if cell_measures_info is None:
             return None
-        
+
         cell_measures_values = {}
 
         for cell_measures_var in cell_measures_info.keys():
 
-            values = deepcopy(cell_measures_info[cell_measures_var])           
+            values = deepcopy(cell_measures_info[cell_measures_var])
             coordinate_len = len(values['data'].shape)
 
             if coordinate_len == 1:
                 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'], 
+                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))
 
             cell_measures_values[cell_measures_var] = values
-            
+
         return cell_measures_values
 
     def _get_lazy_variables(self):
@@ -2101,7 +2107,7 @@ class Nes(object):
 
         # CELL AREA
         if 'cell_area' in self.cell_measures.keys():
-            cell_area = netcdf.createVariable('cell_area', np.float64, self._var_dim, 
+            cell_area = netcdf.createVariable('cell_area', np.float64, self._var_dim,
                                               zlib=self.zip_lvl > 0, complevel=self.zip_lvl)
             if self.size > 1:
                 cell_area.set_collective(True)
@@ -2974,7 +2980,7 @@ class Nes(object):
         x, y, z = pyproj.transform(lla, ecef, lon, lat, np.zeros(lon.shape), radians=False)
 
         return np.column_stack([x, y, z])
-    
+
     def add_4d_vertical_info(self, info_to_add):
         """
         To add the vertical information from other source.
@@ -3011,7 +3017,7 @@ class Nes(object):
             self, new_levels, new_src_vertical=new_src_vertical, kind=kind, extrapolate=extrapolate, info=info)
 
     def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='NearestNeighbour', n_neighbours=4,
-                               info=False, to_providentia=False):
+                               info=False, to_providentia=False, only_create_wm=False, wm=None):
         """
         Horizontal methods from the current grid to another one.
 
@@ -3022,18 +3028,22 @@ class Nes(object):
         weight_matrix_path : str, None
             Path to the weight matrix to read/create.
         kind : str
-            Kind of horizontal methods. choices = ['NearestNeighbour'].
+            Kind of horizontal methods. choices = ['NearestNeighbour', 'Conservative'].
         n_neighbours: int
             Used if kind == NearestNeighbour. Number of nearest neighbours to interpolate. Default: 4.
         info: bool
             Indicates if you want to print extra info during the methods process.
         to_providentia : bool
             Indicates if we want the interpolated grid in Providentia format.
+        only_create_wm : bool
+            Indicates if you want to only create the Weight Matrix.
+        wm : Nes
+            Weight matrix Nes File
         """
         
         return horizontal_interpolation.interpolate_horizontal(
             self, dst_grid, weight_matrix_path=weight_matrix_path, kind=kind, n_neighbours=n_neighbours, info=info, 
-            to_providentia=to_providentia)
+            to_providentia=to_providentia, only_create_wm=only_create_wm, wm=wm)
 
     def calculate_grid_area(self):
         """
@@ -3044,15 +3054,15 @@ class Nes(object):
         self : nes.Nes
             Source projection Nes Object.
         """
-        
+
         grid_area = cell_measures.calculate_grid_area(self)
-        self.cell_measures['cell_area'] = {'data': grid_area.reshape([self.lon_bnds['data'].shape[0], 
+        self.cell_measures['cell_area'] = {'data': grid_area.reshape([self.lon_bnds['data'].shape[0],
                                                                       self.lon_bnds['data'].shape[1]])}
 
         return None
 
     @staticmethod
-    def calculate_geometry_area(geometry_list, earth_radius_minor_axis=6356752.3142, 
+    def calculate_geometry_area(geometry_list, earth_radius_minor_axis=6356752.3142,
                                 earth_radius_major_axis=6378137.0):
         """
         Get coordinate bounds and call function to calculate the area of each cell of a set of geometries.
@@ -3067,7 +3077,7 @@ class Nes(object):
             Radius of the major axis of the Earth.
         """
 
-        return cell_measures.calculate_geometry_area(geometry_list, earth_radius_minor_axis=earth_radius_minor_axis, 
+        return cell_measures.calculate_geometry_area(geometry_list, earth_radius_minor_axis=earth_radius_minor_axis,
                                                      earth_radius_major_axis=earth_radius_major_axis)
 
     @staticmethod
diff --git a/requirements.txt b/requirements.txt
index 3e6c91ae2648806cb55418984623b42ca7a901db..a40856e17a26c8879f41f335d629cc759a16e843 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,8 +1,9 @@
 pycodestyle>=2.10.0
 geopandas>=0.10.2
+rtree>=0.9.0
 pandas>=1.3.5
 netcdf4>=1.6.2
-numpy==1.21.6
+numpy>=1.20.0
 pyproj~=3.2.1
 setuptools>=66.1.1
 pytest>=7.2.1
diff --git a/tests/3.2-test_horiz_interp_bilinear.py b/tests/3.2-test_horiz_interp_bilinear.py
new file mode 100644
index 0000000000000000000000000000000000000000..77bddadbf1daa3fb7655a080068627f534353c3b
--- /dev/null
+++ b/tests/3.2-test_horiz_interp_bilinear.py
@@ -0,0 +1,221 @@
+#!/usr/bin/env python
+import sys
+import timeit
+import pandas as pd
+from mpi4py import MPI
+from nes import *
+import os, sys
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+size = comm.Get_size()
+
+parallel_method = 'T'
+
+result_path = "Times_test_3.2_horiz_interp_bilinear_{0}_{1:03d}.csv".format(parallel_method, size)
+result = pd.DataFrame(index=['read', 'calcul', 'write'],
+                      columns=['3.2.1.Only interp', '3.2.2.Create_WM', "3.2.3.Use_WM", "3.2.4.Read_WM"])
+
+# NAMEE
+src_path = "/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc"
+var_list = ['O3']
+
+# ======================================================================================================================
+# ======================================        Only interp        =====================================================
+# ======================================================================================================================
+test_name = '3.2.1.NN_Only interp'
+if rank == 0:
+    print(test_name)
+
+# READING
+st_time = timeit.default_timer()
+
+# Source data
+src_nes = open_netcdf(src_path, parallel_method=parallel_method)
+src_nes.keep_vars(var_list)
+src_nes.load()
+
+# Destination Grid
+lat_1 = 37
+lat_2 = 43
+lon_0 = -3
+lat_0 = 40
+nx = 397
+ny = 397
+inc_x = 4000
+inc_y = 4000
+x_0 = -807847.688
+y_0 = -797137.125
+dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,
+                     nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method,
+                     times=src_nes.get_full_times())
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+
+interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='NN')
+comm.Barrier()
+result.loc['calcul', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+interp_nes.to_netcdf(test_name.replace(' ', '_') + "{0:03d}.nc".format(size))
+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()
+
+# ======================================================================================================================
+# ======================================         Create_WM        =====================================================
+# ======================================================================================================================
+test_name = '3.2.2.NN_Create_WM'
+if rank == 0:
+    print(test_name)
+
+# READING
+st_time = timeit.default_timer()
+
+# Source data
+src_nes = open_netcdf(src_path, parallel_method=parallel_method)
+src_nes.keep_vars(var_list)
+src_nes.load()
+
+# Destination Grid
+lat_1 = 37
+lat_2 = 43
+lon_0 = -3
+lat_0 = 40
+nx = 397
+ny = 397
+inc_x = 4000
+inc_y = 4000
+x_0 = -807847.688
+y_0 = -797137.125
+dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,
+                     nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method,
+                     times=src_nes.get_full_times())
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+# Cleaning WM
+if os.path.exists("NN_WM_NAMEE_to_IP.nc") and rank == 0:
+    os.remove("NN_WM_NAMEE_to_IP.nc")
+comm.Barrier()
+
+st_time = timeit.default_timer()
+
+wm_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='NN', info=True,
+                                        weight_matrix_path="NN_WM_NAMEE_to_IP.nc", only_create_wm=True)
+comm.Barrier()
+result.loc['calcul', test_name] = timeit.default_timer() - st_time
+
+
+comm.Barrier()
+if rank == 0:
+    print(result.loc[:, test_name])
+sys.stdout.flush()
+
+# ======================================================================================================================
+# ======================================           Use_WM         =====================================================
+# ======================================================================================================================
+test_name = "3.2.3.NN_Use_WM"
+if rank == 0:
+    print(test_name)
+
+# READING
+st_time = timeit.default_timer()
+
+# Source data
+src_nes = open_netcdf(src_path, parallel_method=parallel_method)
+src_nes.keep_vars(var_list)
+src_nes.load()
+
+# Destination Grid
+lat_1 = 37
+lat_2 = 43
+lon_0 = -3
+lat_0 = 40
+nx = 397
+ny = 397
+inc_x = 4000
+inc_y = 4000
+x_0 = -807847.688
+y_0 = -797137.125
+
+dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,
+                     nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method,
+                     times=src_nes.get_full_times())
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+
+interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='NN', wm=wm_nes)
+comm.Barrier()
+result.loc['calcul', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+interp_nes.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size))
+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()
+
+# ======================================================================================================================
+# ======================================           Read_WM         =====================================================
+# ======================================================================================================================
+test_name = "3.2.4.NN_Read_WM"
+if rank == 0:
+    print(test_name)
+
+# READING
+st_time = timeit.default_timer()
+
+# Source data
+src_nes = open_netcdf(src_path, parallel_method=parallel_method)
+src_nes.keep_vars(var_list)
+src_nes.load()
+
+# Destination Grid
+lat_1 = 37
+lat_2 = 43
+lon_0 = -3
+lat_0 = 40
+nx = 397
+ny = 397
+inc_x = 4000
+inc_y = 4000
+x_0 = -807847.688
+y_0 = -797137.125
+
+dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,
+                     nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method,
+                     times=src_nes.get_full_times())
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+
+interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='NN',
+                                            weight_matrix_path="NN_WM_NAMEE_to_IP.nc")
+comm.Barrier()
+result.loc['calcul', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+interp_nes.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size))
+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)
diff --git a/tests/3.3-test_horiz_interp_conservative.py b/tests/3.3-test_horiz_interp_conservative.py
new file mode 100644
index 0000000000000000000000000000000000000000..a6af09955a4690f42548025dd9f2e319805999c4
--- /dev/null
+++ b/tests/3.3-test_horiz_interp_conservative.py
@@ -0,0 +1,227 @@
+#!/usr/bin/env python
+import sys
+import timeit
+import pandas as pd
+from mpi4py import MPI
+from nes import *
+import os, sys
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+size = comm.Get_size()
+
+parallel_method = 'T'
+
+result_path = "Times_test_3.3_horiz_interp_conservative.py_{0}_{1:03d}.csv".format(parallel_method, size)
+result = pd.DataFrame(index=['read', 'calcul', 'write'],
+                      columns=['3.3.1.Only interp', '3.3.2.Create_WM', "3.3.3.Use_WM", "3.3.4.Read_WM"])
+
+# NAMEE
+src_path = "/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc"
+var_list = ['O3']
+
+# ======================================================================================================================
+# ======================================        Only interp        =====================================================
+# ======================================================================================================================
+test_name = '3.3.1.Only interp'
+if rank == 0:
+    print(test_name)
+
+# READING    final_dst.variables[var_name]['data'][time, lev] = np.sum(weights * src_aux, axis=1)
+
+st_time = timeit.default_timer()
+
+# Source data
+src_nes = open_netcdf(src_path, parallel_method=parallel_method)
+src_nes.keep_vars(var_list)
+src_nes.load()
+
+# Destination Grid
+lat_1 = 37
+lat_2 = 43
+lon_0 = -3
+lat_0 = 40
+nx = 397
+ny = 397
+inc_x = 4000
+inc_y = 4000
+x_0 = -807847.688
+y_0 = -797137.125
+dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,
+                     nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method,
+                     times=src_nes.get_full_times())
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+
+interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative')
+# interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', weight_matrix_path='T_WM.nc')
+
+comm.Barrier()
+result.loc['calcul', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+interp_nes.to_netcdf(test_name.replace(' ', '_') + "{0:03d}.nc".format(size), serial=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()
+
+# ======================================================================================================================
+# ======================================         Create_WM        =====================================================
+# ======================================================================================================================
+test_name = '3.3.2.Create_WM'
+if rank == 0:
+    print(test_name)
+
+# READING
+st_time = timeit.default_timer()
+
+# Source data
+src_nes = open_netcdf(src_path, parallel_method=parallel_method)
+src_nes.keep_vars(var_list)
+src_nes.load()
+
+# Destination Grid
+lat_1 = 37
+lat_2 = 43
+lon_0 = -3
+lat_0 = 40
+nx = 397
+ny = 397
+inc_x = 4000
+inc_y = 4000
+x_0 = -807847.688
+y_0 = -797137.125
+dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,
+                     nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method,
+                     times=src_nes.get_full_times())
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+# Cleaning WM
+if os.path.exists("CONS_WM_NAMEE_to_IP.nc") and rank == 0:
+    os.remove("CONS_WM_NAMEE_to_IP.nc")
+comm.Barrier()
+
+st_time = timeit.default_timer()
+
+wm_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', info=True,
+                                        weight_matrix_path="CONS_WM_NAMEE_to_IP.nc", only_create_wm=True)
+comm.Barrier()
+result.loc['calcul', test_name] = timeit.default_timer() - st_time
+
+# st_time = timeit.default_timer()
+# interp_nes.to_netcdf(test_name.replace(' ', '_') + ".nc")
+# 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()
+
+# ======================================================================================================================
+# ======================================           Use_WM         =====================================================
+# ======================================================================================================================
+test_name = "3.3.3.Use_WM"
+if rank == 0:
+    print(test_name)
+
+# READING
+st_time = timeit.default_timer()
+
+# Source data
+src_nes = open_netcdf(src_path, parallel_method=parallel_method)
+src_nes.keep_vars(var_list)
+src_nes.load()
+
+# Destination Grid
+lat_1 = 37
+lat_2 = 43
+lon_0 = -3
+lat_0 = 40
+nx = 397
+ny = 397
+inc_x = 4000
+inc_y = 4000
+x_0 = -807847.688
+y_0 = -797137.125
+
+dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,
+                     nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method,
+                     times=src_nes.get_full_times())
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+
+interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', wm=wm_nes)
+comm.Barrier()
+result.loc['calcul', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+interp_nes.to_netcdf(test_name.replace(' ', '_') + "{0:03d}.nc".format(size))
+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()
+
+# ======================================================================================================================
+# ======================================           Read_WM         =====================================================
+# ======================================================================================================================
+test_name = "3.3.4.Read_WM"
+if rank == 0:
+    print(test_name)
+
+# READING
+st_time = timeit.default_timer()
+
+# Source data
+src_nes = open_netcdf(src_path, parallel_method=parallel_method)
+src_nes.keep_vars(var_list)
+src_nes.load()
+
+# Destination Grid
+lat_1 = 37
+lat_2 = 43
+lon_0 = -3
+lat_0 = 40
+nx = 397
+ny = 397
+inc_x = 4000
+inc_y = 4000
+x_0 = -807847.688
+y_0 = -797137.125
+
+dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,
+                     nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, parallel_method=parallel_method,
+                     times=src_nes.get_full_times())
+comm.Barrier()
+result.loc['read', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+
+interp_nes = src_nes.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', weight_matrix_path="CONS_WM_NAMEE_to_IP.nc")
+comm.Barrier()
+result.loc['calcul', test_name] = timeit.default_timer() - st_time
+
+st_time = timeit.default_timer()
+interp_nes.to_netcdf(test_name.replace(' ', '_') + "{0:03d}.nc".format(size))
+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)
diff --git a/tests/sbatch_test_MN4.cmd b/tests/sbatch_test_MN4.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..09011158b106e8716aab17388ac4c3e5a26c1569
--- /dev/null
+++ b/tests/sbatch_test_MN4.cmd
@@ -0,0 +1,35 @@
+#!/bin/bash
+
+#SBATCH --qos=debug
+#SBATCH -A bsc32
+#SBATCH --cpus-per-task=1
+#SBATCH -n 1
+#SBATCH -t 00:30:00
+#SBATCH -J NES_mn4_tests
+#SBATCH --output=log_mn4_NES_%j.out
+#SBATCH --error=log_mn4_NES_%j.err
+#SBATCH --exclusive
+
+### ulimit -s 128000
+
+module purge
+module use /gpfs/projects/bsc32/software/suselinux/11/modules/all
+
+module load Python/3.7.4-GCCcore-8.3.0
+module load netcdf4-python/1.5.3-foss-2019b-Python-3.7.4
+module load cfunits/1.8-foss-2019b-Python-3.7.4
+module load xarray/0.17.0-foss-2019b-Python-3.7.4
+module load OpenMPI/4.0.5-GCC-8.3.0-mn4
+module load filelock/3.7.1-foss-2019b-Python-3.7.4
+module load eccodes-python/0.9.5-foss-2019b-Python-3.7.4
+module load pyproj/2.5.0-foss-2019b-Python-3.7.4
+module load geopandas/0.8.1-foss-2019b-Python-3.7.4
+module load Shapely/1.7.1-foss-2019b-Python-3.7.4
+
+# export PYTHONPATH=/gpfs/projects/bsc32/models/NES_master:${PYTHONPATH}
+# cd /gpfs/projects/bsc32/models/NES_master/tests
+
+export PYTHONPATH=/gpfs/scratch/bsc32/bsc32538/NES_tests/NES:${PYTHONPATH}
+cd /gpfs/scratch/bsc32/bsc32538/NES_tests/NES/tests
+
+mpirun --mca mpi_warn_on_fork 0 -np 1 python test_interp_conservative.py
diff --git a/tests/test_bash_mn4.cmd b/tests/test_bash_mn4.cmd
index 5a0fca92696b2fd4329d479b3f122dc5e98d700d..9dfd27dbb3826faa0805dc9d1a552a5fcf81424f 100644
--- a/tests/test_bash_mn4.cmd
+++ b/tests/test_bash_mn4.cmd
@@ -24,3 +24,4 @@ cd /gpfs/projects/bsc32/models/NES_master/tests
 #mpirun --mca mpi_warn_on_fork 0 -np 4 python 2.1-test_spatial_join.py
 #mpirun --mca mpi_warn_on_fork 0 -np 4 python 2.3-test_bounds.py
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 2.4-test_cell_area.py
+mpirun --mca mpi_warn_on_fork 0 -np 4 python 5.3-test_horiz_interp_conservative.py
diff --git a/tests/test_bash_nord3v2.cmd b/tests/test_bash_nord3v2.cmd
index 26d9e28448284eb5191787a3031f1627dc446f67..e585eab22d231734eae8dc1545e973ab200c8ae0 100644
--- a/tests/test_bash_nord3v2.cmd
+++ b/tests/test_bash_nord3v2.cmd
@@ -23,3 +23,4 @@ cd /gpfs/projects/bsc32/models/NES_master/tests
 #mpirun --mca mpi_warn_on_fork 0 -np 4 python 2.1-test_spatial_join.py
 #mpirun --mca mpi_warn_on_fork 0 -np 4 python 2.3-test_bounds.py
 mpirun --mca mpi_warn_on_fork 0 -np 4 python 2.4-test_cell_area.py
+mpirun --mca mpi_warn_on_fork 0 -np 4 python 5.3-test_horiz_interp_conservative.py
diff --git a/tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb b/tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb
index 38b1ac5812fc05899337fd3a12ed2bca05023ee9..d126994d0a6d8a0b44c709d6838a1a8954fdc436 100644
--- a/tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb
+++ b/tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# How to interpolate vertically"
+    "# Vertical interpolation"
    ]
   },
   {
diff --git a/tutorials/4.Interpolation/4.2.Horizontal_Interpolation.ipynb b/tutorials/4.Interpolation/4.2.Horizontal_Interpolation.ipynb
index 5f0379f37316bb8f5f746651e8218948aa248f6d..63264addac257ca827dd7f52091318bc36e8399c 100644
--- a/tutorials/4.Interpolation/4.2.Horizontal_Interpolation.ipynb
+++ b/tutorials/4.Interpolation/4.2.Horizontal_Interpolation.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# How to interpolate horizontally"
+    "# Horizontal interpolation"
    ]
   },
   {
@@ -593,7 +593,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8 (default, Aug 12 2021, 07:06:15) \n[GCC 8.4.1 20200928 (Red Hat 8.4.1-1)]"
+   "version": "3.7.4"
   },
   "vscode": {
    "interpreter": {
diff --git a/tutorials/4.Interpolation/4.3.Conservative_Interpolation.ipynb b/tutorials/4.Interpolation/4.3.Conservative_Interpolation.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..6177f4cc964a8b92a7fabb9eeeea7a15acb466fc
--- /dev/null
+++ b/tutorials/4.Interpolation/4.3.Conservative_Interpolation.ipynb
@@ -0,0 +1,396 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Horizontal conservative interpolation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from nes import *"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 1. Read data to interpolate"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Original path: /esarchive/exp/monarch/a4dd/original_files/000/2022111512/MONARCH_d01_2022111512.nc\n",
+    "# Rotated grid from MONARCH (NAMEE)\n",
+    "source_path = \"/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "source_grid = open_netcdf(path=source_path, info=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'data': masked_array(\n",
+       "   data=[[16.350338, 16.43293 , 16.515146, ..., 16.515146, 16.43293 ,\n",
+       "          16.350338],\n",
+       "         [16.527426, 16.610239, 16.692677, ..., 16.692677, 16.610243,\n",
+       "          16.527426],\n",
+       "         [16.704472, 16.787508, 16.870167, ..., 16.870167, 16.78751 ,\n",
+       "          16.704472],\n",
+       "         ...,\n",
+       "         [58.32095 , 58.472683, 58.62431 , ..., 58.62431 , 58.472683,\n",
+       "          58.32095 ],\n",
+       "         [58.426285, 58.5782  , 58.730026, ..., 58.730026, 58.5782  ,\n",
+       "          58.426285],\n",
+       "         [58.530792, 58.6829  , 58.83492 , ..., 58.83492 , 58.682903,\n",
+       "          58.530792]],\n",
+       "   mask=False,\n",
+       "   fill_value=1e+20,\n",
+       "   dtype=float32),\n",
+       " 'dimensions': ('rlat', 'rlon'),\n",
+       " 'long_name': 'latitude',\n",
+       " 'units': 'degrees_north',\n",
+       " 'standard_name': 'latitude',\n",
+       " 'coordinates': 'lon lat'}"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "source_grid.lat"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'data': masked_array(\n",
+       "   data=[[-22.181265, -22.016672, -21.851799, ...,  41.851795,  42.016666,\n",
+       "           42.18126 ],\n",
+       "         [-22.27818 , -22.113186, -21.947905, ...,  41.9479  ,  42.113174,\n",
+       "           42.27817 ],\n",
+       "         [-22.375267, -22.209873, -22.04419 , ...,  42.044186,  42.209873,\n",
+       "           42.375263],\n",
+       "         ...,\n",
+       "         [-67.57767 , -67.397064, -67.21535 , ...,  87.21534 ,  87.39706 ,\n",
+       "           87.57766 ],\n",
+       "         [-67.90188 , -67.72247 , -67.54194 , ...,  87.54194 ,  87.72246 ,\n",
+       "           87.90187 ],\n",
+       "         [-68.228035, -68.04982 , -67.870514, ...,  87.87051 ,  88.04982 ,\n",
+       "           88.228035]],\n",
+       "   mask=False,\n",
+       "   fill_value=1e+20,\n",
+       "   dtype=float32),\n",
+       " 'dimensions': ('rlat', 'rlon'),\n",
+       " 'long_name': 'longitude',\n",
+       " 'units': 'degrees_east',\n",
+       " 'standard_name': 'longitude',\n",
+       " 'coordinates': 'lon lat'}"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "source_grid.lon"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Rank 000: Loading O3 var (1/1)\n",
+      "Rank 000: Loaded O3 var ((37, 24, 271, 351))\n"
+     ]
+    }
+   ],
+   "source": [
+    "source_grid.keep_vars('O3')\n",
+    "source_grid.load()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 2. Create destination grid"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lat_1 = 37\n",
+    "lat_2 = 43\n",
+    "lon_0 = -3\n",
+    "lat_0 = 40\n",
+    "nx = 397\n",
+    "ny = 397\n",
+    "inc_x = 4000\n",
+    "inc_y = 4000\n",
+    "x_0 = -807847.688\n",
+    "y_0 = -797137.125\n",
+    "dst_nes = create_nes(comm=None, info=False, projection='lcc', lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0,\n",
+    "                     nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0, times=source_grid.get_full_times())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 3. Interpolation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 3.1. Without creating weight matrix"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "interp_nes = source_grid.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'data': array([[32.49460816, 32.49803615, 32.50144726, ..., 32.52460207,\n",
+       "         32.52130788, 32.51799681],\n",
+       "        [32.53024662, 32.5336765 , 32.5370895 , ..., 32.5602571 ,\n",
+       "         32.55696109, 32.55364819],\n",
+       "        [32.5658877 , 32.56931947, 32.57273435, ..., 32.59591474,\n",
+       "         32.59261692, 32.58930218],\n",
+       "        ...,\n",
+       "        [46.60231473, 46.60653579, 46.6107361 , ..., 46.63924881,\n",
+       "         46.63519229, 46.63111499],\n",
+       "        [46.63791957, 46.64214277, 46.6463452 , ..., 46.67487234,\n",
+       "         46.67081377, 46.66673441],\n",
+       "        [46.67352141, 46.67774674, 46.68195131, ..., 46.71049289,\n",
+       "         46.70643226, 46.70235083]])}"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "interp_nes.lat"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'data': array([[-11.56484687, -11.52259289, -11.48033507, ...,   5.18764453,\n",
+       "           5.22992847,   5.27220869],\n",
+       "        [-11.56892309, -11.52664925, -11.48437156, ...,   5.1915433 ,\n",
+       "           5.23384715,   5.27614727],\n",
+       "        [-11.57300318, -11.53070946, -11.48841188, ...,   5.19544578,\n",
+       "           5.23776955,   5.2800896 ],\n",
+       "        ...,\n",
+       "        [-13.53865445, -13.48682654, -13.4349915 , ...,   7.07590089,\n",
+       "           7.12778439,   7.17966098],\n",
+       "        [-13.54481681, -13.49295916, -13.44109437, ...,   7.08179741,\n",
+       "           7.13371074,   7.18561716],\n",
+       "        [-13.55098635, -13.49909893, -13.44720434, ...,   7.0877008 ,\n",
+       "           7.139644  ,   7.19158028]])}"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "interp_nes.lon"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 3.2. Creating weight matrix"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Creating Weight Matrix\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/esarchive/scratch/avilanova/software/NES/nes/nc_projections/default_nes.py:2035: DeprecationWarning: tostring() is deprecated. Use tobytes() instead.\n",
+      "  time_var.units, time_var.calendar)\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Weight Matrix done!\n"
+     ]
+    }
+   ],
+   "source": [
+    "wm_nes = source_grid.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', info=True,\n",
+    "                                            weight_matrix_path=\"CONS_WM_NAMEE_to_IP.nc\", only_create_wm=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "interp_nes = source_grid.interpolate_horizontal(dst_grid=dst_nes, kind='Conservative', wm=wm_nes)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'data': array([[32.49460816, 32.49803615, 32.50144726, ..., 32.52460207,\n",
+       "         32.52130788, 32.51799681],\n",
+       "        [32.53024662, 32.5336765 , 32.5370895 , ..., 32.5602571 ,\n",
+       "         32.55696109, 32.55364819],\n",
+       "        [32.5658877 , 32.56931947, 32.57273435, ..., 32.59591474,\n",
+       "         32.59261692, 32.58930218],\n",
+       "        ...,\n",
+       "        [46.60231473, 46.60653579, 46.6107361 , ..., 46.63924881,\n",
+       "         46.63519229, 46.63111499],\n",
+       "        [46.63791957, 46.64214277, 46.6463452 , ..., 46.67487234,\n",
+       "         46.67081377, 46.66673441],\n",
+       "        [46.67352141, 46.67774674, 46.68195131, ..., 46.71049289,\n",
+       "         46.70643226, 46.70235083]])}"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "interp_nes.lat"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'data': array([[-11.56484687, -11.52259289, -11.48033507, ...,   5.18764453,\n",
+       "           5.22992847,   5.27220869],\n",
+       "        [-11.56892309, -11.52664925, -11.48437156, ...,   5.1915433 ,\n",
+       "           5.23384715,   5.27614727],\n",
+       "        [-11.57300318, -11.53070946, -11.48841188, ...,   5.19544578,\n",
+       "           5.23776955,   5.2800896 ],\n",
+       "        ...,\n",
+       "        [-13.53865445, -13.48682654, -13.4349915 , ...,   7.07590089,\n",
+       "           7.12778439,   7.17966098],\n",
+       "        [-13.54481681, -13.49295916, -13.44109437, ...,   7.08179741,\n",
+       "           7.13371074,   7.18561716],\n",
+       "        [-13.55098635, -13.49909893, -13.44720434, ...,   7.0877008 ,\n",
+       "           7.139644  ,   7.19158028]])}"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "interp_nes.lon"
+   ]
+  }
+ ],
+ "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
+}
diff --git a/tutorials/4.Interpolation/4.3.Providentia_Interpolation.ipynb b/tutorials/4.Interpolation/4.4.Providentia_Interpolation.ipynb
similarity index 100%
rename from tutorials/4.Interpolation/4.3.Providentia_Interpolation.ipynb
rename to tutorials/4.Interpolation/4.4.Providentia_Interpolation.ipynb
diff --git a/tutorials/4.Interpolation/4.4.NES_vs_Providentia_Interpolation.ipynb b/tutorials/4.Interpolation/4.5.NES_vs_Providentia_Interpolation.ipynb
similarity index 100%
rename from tutorials/4.Interpolation/4.4.NES_vs_Providentia_Interpolation.ipynb
rename to tutorials/4.Interpolation/4.5.NES_vs_Providentia_Interpolation.ipynb