diff --git a/inputs/namelist_cmems.bdy b/inputs/namelist_cmems.bdy
index d57ecce8ea03e8bcfd64c4cfc5a43eb144951ce7..2bc9c3aa8818a87126818e53b3019eeb1c742eea 100644
--- a/inputs/namelist_cmems.bdy
+++ b/inputs/namelist_cmems.bdy
@@ -52,7 +52,7 @@
 !  CMEMS Data Source Configuration
 !------------------------------------------------------------------------------
    ln_use_cmems             = .true.
-   ln_download_cmems        = .true.
+   ln_download_cmems        = .false.
    sn_cmems_dir             = '/Users/thopri/Projects/PyNEMO/inputs/' ! where to download CMEMS input files (static and variable)
    ln_download_static       = .false.
    ln_subset_static         = .false.
@@ -103,7 +103,7 @@
 !------------------------------------------------------------------------------
 !  unstructured open boundaries tidal parameters                        
 !------------------------------------------------------------------------------
-    ln_tide        = .false.              !  =T : produce bdy tidal conditions
+    ln_tide        = .true.              !  =T : produce bdy tidal conditions
     sn_tide_model  = 'fes'                !  Name of tidal model (fes|tpxo)
     clname(1)      = 'M2'                 !  constituent name
     clname(2)      = 'S2'
diff --git a/pynemo/tide/fes_extract_HC.py b/pynemo/tide/fes_extract_HC.py
index 56508558d3bcef0c73aecca64362c34c2568c3ae..832bdbaab3aefb2d7c272331e3cc03d4cb9afc4c 100644
--- a/pynemo/tide/fes_extract_HC.py
+++ b/pynemo/tide/fes_extract_HC.py
@@ -53,26 +53,20 @@ class HcExtract(object):
             lon_resolution = lon_z[1] - lon_z[0]
             data_in_km = 0 # added to maintain the reference to matlab tmd code
 
-            # extract example amplitude grid for Z, U and V and change NaNs to 0 (for land) and other values to 1 (for water)
-            mask_z = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['amplitude'][:]))
-            where_are_NaNs = np.isnan(mask_z)
-            where_no_NaNs = np.invert(np.isnan(mask_z))
-            mask_z[where_no_NaNs] = 1
-            mask_z[where_are_NaNs] = 0
+            # extract example amplitude grid for Z, U and V and change NaNs to 1 (for land) and other values to 0 (for water)
+            mask_z = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['amplitude'][:])))
+            mask_z[mask_z == 18446744073709551616.00000] = 0
+            mask_z[mask_z != 0] = 1
             self.mask_dataset[mz_name] = mask_z
 
-            mask_u = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_U.nc').variables['Ua'][:]))
-            where_are_NaNs = np.isnan(mask_u)
-            where_no_NaNs = np.invert(np.isnan(mask_u))
-            mask_u[where_no_NaNs] = 1
-            mask_u[where_are_NaNs] = 0
+            mask_u = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_U.nc').variables['Ua'][:])))
+            mask_u[mask_u == 18446744073709551616.00000] = 0
+            mask_u[mask_u != 0] = 1
             self.mask_dataset[mu_name] = mask_u
 
-            mask_v = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_V.nc').variables['Va'][:]))
-            where_are_NaNs = np.isnan(mask_v)
-            where_no_NaNs = np.invert(np.isnan(mask_v))
-            mask_v[where_no_NaNs] = 1
-            mask_v[where_are_NaNs] = 0
+            mask_v = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_V.nc').variables['Va'][:])))
+            mask_v[mask_v == 18446744073709551616.00000] = 0
+            mask_v[mask_v != 0] = 1
             self.mask_dataset[mv_name] = mask_v
 
 
@@ -82,8 +76,10 @@ class HcExtract(object):
             lat_z = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['lat'][:])
             lon_z = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['lon'][:])
             for ncon in range(len(constituents)):
-                amp = np.array(np.rot90(Dataset(settings['tide_fes']+str(constituents[ncon])+'_Z.nc').variables['amplitude'][:]))
-                phase = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:]))
+                amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+str(constituents[ncon])+'_Z.nc').variables['amplitude'][:])))
+                amp[amp == 18446744073709551616.00000] = 0
+                amp = amp/100.00
+                phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:])))
                 hRe.append(amp*np.sin(phase))
                 hIm.append(amp*np.cos(phase))
             hRe = np.stack(hRe)
@@ -96,8 +92,10 @@ class HcExtract(object):
             lat_u = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['lat'][:])
             lon_u = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['lon'][:])
             for ncon in range(len(constituents)):
-                amp = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:]))
-                phase = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:]))
+                amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:])))
+                amp[amp == 18446744073709551616.00000] = 0
+                amp = amp/100.00
+                phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:])))
                 URe.append(amp*np.sin(phase))
                 UIm.append(amp*np.cos(phase))
             URe = np.stack(URe)
@@ -109,8 +107,10 @@ class HcExtract(object):
             lat_v = np.array(Dataset(settings['tide_fes'] + constituents[ncon] + '_V.nc').variables['lat'][:])
             lon_v = np.array(Dataset(settings['tide_fes'] + constituents[ncon] + '_V.nc').variables['lon'][:])
             for ncon in range(len(constituents)):
-                amp = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:]))
-                phase = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:]))
+                amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:])))
+                amp[amp == 18446744073709551616.00000] = 0
+                amp = amp/100.00
+                phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:])))
                 VRe.append(amp*np.sin(phase))
                 VIm.append(amp*np.cos(phase))
             VRe = np.stack(VRe)
diff --git a/test_scripts/bdy_var_plot.py b/test_scripts/bdy_var_plot.py
index be1e8136ca3fe4fc2f585a20e90b59b8f6c4ebcd..0d8bae13693ed78b23ff98bd9a7966f0f71da249 100644
--- a/test_scripts/bdy_var_plot.py
+++ b/test_scripts/bdy_var_plot.py
@@ -61,13 +61,13 @@ def nemo_bdy_order(fname):
         id_order[0,] = 0
         flag = False
         mark = 0
-        source_tree = sp.cKDTree(zip(nbi_tmp, nbj_tmp), balanced_tree=False, compact_nodes=False)
+        source_tree = sp.cKDTree(list(zip(nbi_tmp, nbj_tmp)), balanced_tree=False, compact_nodes=False)
 
         # order bdy entries
 
         while count < nbdy[r]:
 
-            nn_dist, nn_id = source_tree.query(zip(nbi_tmp[id_order[count - 1]], nbj_tmp[id_order[count - 1]]),
+            nn_dist, nn_id = source_tree.query(list(zip(nbi_tmp[id_order[count - 1]], nbj_tmp[id_order[count - 1]])),
                                                k=3, distance_upper_bound=2.9)
             if np.sum(id_order == nn_id[0, 1]) == 1:  # is the nearest point already in the list?
                 if np.sum(id_order == nn_id[0, 2]) == 1:  # is the 2nd nearest point already in the list?
@@ -207,19 +207,19 @@ def plot_bdy(fname, bdy_ind, bdy_dst, bdy_brk, varnam, t, rw):
 
         # create a pseudo bathymetry from the depth data
 
-        bathy = np.zeros_like(coords)
-        mbath = np.sum(dta[n].mask == 0, axis=0)
+        #bathy = np.zeros_like(coords)
+        #mbath = np.sum(dta[n].mask == 0, axis=0)
 
-        for i in range(len(coords)):
-            bathy[i] = gdepw[mbath[i], i]
+        #for i in range(len(coords)):
+        #    bathy[i] = gdepw[mbath[i], i]
 
-        bathy_patch = Polygon(np.vstack((np.hstack((coords[0], coords, coords[-1])),
-                                         np.hstack((np.amax(bathy[:]), bathy, np.amax(bathy[:]))))).T,
-                              closed=True,
-                              facecolor=(0.8, 0.8, 0), alpha=0, edgecolor=None)
+        #bathy_patch = Polygon(np.vstack((np.hstack((coords[0], coords, coords[-1])),
+        #                                 np.hstack((np.amax(bathy[:]), bathy, np.amax(bathy[:]))))).T,
+        #                      closed=True,
+        #                      facecolor=(0.8, 0.8, 0), alpha=0, edgecolor=None)
 
         # Add patch to axes
-        ax[n].add_patch(bathy_patch)
+        #ax[n].add_patch(bathy_patch)
         ax[n].set_title('BDY points along section: ' + str(n))
         patches = []
         colors = []
@@ -248,17 +248,17 @@ def plot_bdy(fname, bdy_ind, bdy_dst, bdy_brk, varnam, t, rw):
         #         plt.plot(x, y, 'k-', linewidth=0.1)
         #         plt.plot(coords[i], gdept[k, i], 'k.', markersize=1)
 
-        plt.plot(coords, bathy, '-', color=(0.4, 0, 0))
+        #plt.plot(coords, bathy, '-', color=(0.4, 0, 0))
         p = PatchCollection(patches, alpha=0.8, edgecolor='none')
         p.set_array(np.array(colors))
         ax[n].add_collection(p)
         f.colorbar(p, ax=ax[n])
-        ax[n].set_ylim((0, np.max(bathy)))
+        #ax[n].set_ylim((0, np.max(bathy)))
         ax[n].invert_yaxis()
 
     return f
 
-fname = '/Users/thopri/Projects/PyNEMO/outputs/NNA_R12_bdyT_y2017m11.nc'
+fname = '/Users/thopri/Projects/PyNEMO/outputs/NNA_R12_bdyT_y2017m01.nc'
 print(fname)
 ind, dst, brk = nemo_bdy_order(fname)
 f = plot_bdy(fname, ind, dst, brk, 'thetao', 0, 0)
diff --git a/test_scripts/test_cmems.py b/test_scripts/test_cmems.py
deleted file mode 100644
index 75362698206ddca2e0a665f6f6b16ea425c3ea53..0000000000000000000000000000000000000000
--- a/test_scripts/test_cmems.py
+++ /dev/null
@@ -1,97 +0,0 @@
-#!/usr/bin/env python2
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Oct 31 15:49:52 2019
-
-Test script to test download CMEMS data function in PyNEMO. To be incorprated
-into profile.py
-
-@author: thopri
-"""
-from pynemo import nemo_bdy_dl_cmems as dl_cmems
-from pynemo import nemo_bdy_setup as setup
-from calendar import monthrange
-import sys
-import time
-import threading
-
-Setup = setup.Setup('/Users/thopri/Projects/PyNEMO/inputs/namelist_cmems.bdy') # default settings file
-settings = Setup.settings
-
-if settings['use_cmems'] == True:
-
-    if settings['year_end'] - settings['year_000'] > 0:
-        date_min = settings['year_000']+'-01-01'
-        date_max = settings['year_end']+'-12-31'
-
-    days_mth = monthrange(settings['year_end'],settings['month_end'])
-
-    date_min = str(settings['year_000'])+'-'+str(settings['month_000'])+'-01'
-
-    date_max = str(settings['year_end'])+'-'+str(settings['month_end'])+'-'+str(days_mth[1])
-
-    cmes_config= {
-           'ini_config_template'    : settings['cmes_config_template'],
-           'user'                   : settings['cmes_usr'],
-           'pwd'                    : settings['cmes_pwd'],
-           'motu_server'            : settings['motu_server'],
-           'service_id'             : settings['cmes_model'],
-           'product_id'             : settings['cmes_product'],
-           'date_min'               : date_min,
-           'date_max'               : date_max,
-           'latitude_min'           : settings['latitude_min'],
-           'latitude_max'           : settings['latitude_max'],
-           'longitude_min'          : settings['longitude_min'],
-           'longitude_max'          : settings['longitude_max'],
-           'depth_min'              : settings['depth_min'],
-           'depth_max'              : settings['depth_max'],
-           'variable'               : settings['cmes_variable'],
-           'src_dir'                : settings['src_dir'],
-           'out_name'               : settings['cmes_output'],
-           'config_out'             : settings['cmes_config']
-           }
-
-    chk = dl_cmems.chk_motu()
-
-    dl = dl_cmems.request_CMEMS(cmes_config)
-
-class progress_bar_loading(threading.Thread):
-
-    def run(self):
-            global stop
-            global kill
-            print 'Loading....  ',
-            sys.stdout.flush()
-            i = 0
-            while stop != True:
-                    if (i%4) == 0:
-                        sys.stdout.write('\b/')
-                    elif (i%4) == 1:
-                        sys.stdout.write('\b-')
-                    elif (i%4) == 2:
-                        sys.stdout.write('\b\\')
-                    elif (i%4) == 3:
-                        sys.stdout.write('\b|')
-
-                    sys.stdout.flush()
-                    time.sleep(0.2)
-                    i+=1
-
-            if kill == True:
-                print '\b\b\b\b ABORT!',
-            else:
-                print '\b\b done!',
-
-
-kill = False
-stop = False
-p = progress_bar_loading()
-p.start()
-
-try:
-    #anything you want to run.
-    time.sleep(1)
-    stop = True
-except KeyboardInterrupt or EOFError:
-         kill = True
-         stop = True
\ No newline at end of file