diff --git a/inputs/namelist_remote.bdy b/inputs/namelist_remote.bdy
index 72b58c848840f7985bcd296a2fc4f6b475442701..bb2a1a78e12be10461347cdd940b3bb3483b6bb1 100644
--- a/inputs/namelist_remote.bdy
+++ b/inputs/namelist_remote.bdy
@@ -56,7 +56,7 @@
     ln_mask_file   = .false.              !  =T : read mask from file
     cn_mask_file   = 'mask.nc'            !  name of mask file 
                                           !  (if ln_mask_file=.TRUE.)
-    ln_dyn2d       = .false.              !  boundary conditions for 
+    ln_dyn2d       = .true.              !  boundary conditions for
                                           !  barotropic fields
     ln_dyn3d       = .false.              !  boundary conditions for 
                                           !  baroclinic velocities
diff --git a/pynemo/nemo_bdy_extr_tm3.py b/pynemo/nemo_bdy_extr_tm3.py
index 632e03573fc4dffafc87b7d48e0a95e0ae28dc61..584ff48b46d82e139d385976f25d8b3e3670a1b7 100644
--- a/pynemo/nemo_bdy_extr_tm3.py
+++ b/pynemo/nemo_bdy_extr_tm3.py
@@ -894,12 +894,12 @@ class Extract:
             
 #        for v in self.variables:
         for v in self.var_nam:
-            if self.settings['dyn2d']: # Calculate depth averaged velocity
-                tile_dz = np.tile(self.bdy_dz, [len(self.time_counter), 1, 1, 1])
-                tmp_var = np.reshape(self.d_bdy[v][year]['data'][:,:,:], tile_dz.shape)
-                tmp_var = np.nansum(tmp_var * tile_dz, 2) /np.nansum(tile_dz, 2)
-            else: # Replace NaNs with specified fill value
-                tmp_var = np.where(np.isnan(self.d_bdy[v][year]['data'][:,:,:]),
+            #if self.settings['dyn2d']: # Calculate depth averaged velocity
+            #    tile_dz = np.tile(self.bdy_dz, [len(self.time_counter), 1, 1, 1])
+            #    tmp_var = np.reshape(self.d_bdy[v][year]['data'][:,:,:], tile_dz.shape)
+            #    tmp_var = np.nansum(tmp_var * tile_dz, 2) /np.nansum(tile_dz, 2)
+            #else: # Replace NaNs with specified fill value
+            tmp_var = np.where(np.isnan(self.d_bdy[v][year]['data'][:,:,:]),
                                             self.settings['fv'], 
                                             self.d_bdy[v][year]['data'][:,:,:])
                
diff --git a/pynemo/nemo_bdy_ncgen.py b/pynemo/nemo_bdy_ncgen.py
index 96e1e95c9c66c43e8bb24e34b4ded005ade37d69..c0ccccc2e3de2bacf0c75c05fe8afe2b61ebf738 100644
--- a/pynemo/nemo_bdy_ncgen.py
+++ b/pynemo/nemo_bdy_ncgen.py
@@ -240,7 +240,7 @@ def CreateBDYNetcdfFile(filename, N, I, J, K, rw, h, orig, fv, calendar, grd, va
 
         ncid.close()
 
-    if var_nam[0] == 'thetao' or var_nam[0] == 'so':
+    if var_nam[0] == 'thetao' or var_nam[0] == 'so' or var_nam[0] == 'uo' or var_nam[0] == 'vo':
         logging.info('CMEMS variables identified, using default CMEMS variables.....')
         gridNames = ['T', 'I', 'U', 'V', 'E', 'Z'] # All possible grids
 
@@ -286,6 +286,7 @@ def CreateBDYNetcdfFile(filename, N, I, J, K, rw, h, orig, fv, calendar, grd, va
             vartmpID = ncid.createVariable('thetao', 'f4', ('time_counter', 'z', 'yb', 'xb', ),
                                           fill_value=fv)
             varsalID = ncid.createVariable('so', 'f4', ('time_counter', 'z', 'yb', 'xb', ), fill_value=fv)
+            varsshID = ncid.createVariable('zos', 'f4', ('time_counter', 'z', 'yb', 'xb',), fill_value=fv)
             if grd == 'I':
                 varildID = ncid.createVariable('ileadfra', 'f4', ('time_counter', 'yb', 'xb',),
                                                fill_value=fv)
@@ -297,16 +298,16 @@ def CreateBDYNetcdfFile(filename, N, I, J, K, rw, h, orig, fv, calendar, grd, va
             varztID = ncid.createVariable('depthu', 'f4', ('z', 'yb', 'xb', ), fill_value=fv)
             varbtuID = ncid.createVariable('vobtcrtx', 'f4', ('time_counter', 'yb', 'xb', ),
                                            fill_value=fv)
-            vartouID = ncid.createVariable('vozocrtx', 'f4', ('time_counter', 'z', 'yb', 'xb', ),
+            vartouID = ncid.createVariable('uo', 'f4', ('time_counter', 'z', 'yb', 'xb', ),
                                            fill_value=fv)
         elif grd == 'V':
             varztID = ncid.createVariable('depthv', 'f4', ('z', 'yb', 'xb', ))
             varbtvID = ncid.createVariable('vobtcrty', 'f4', ('time_counter', 'yb', 'xb', ),
                                            fill_value=fv)
-            vartovID = ncid.createVariable('vomecrty', 'f4', ('time_counter', 'z', 'yb', 'xb',),
+            vartovID = ncid.createVariable('vo', 'f4', ('time_counter', 'z', 'yb', 'xb',),
                                            fill_value=fv)
         elif grd == 'Z':
-            varsshID = ncid.createVariable('sossheig', 'f4', ('time_counter', 'yb', 'xb', ),
+            varsshID = ncid.createVariable('zos', 'f4', ('time_counter', 'yb', 'xb', ),
                                            fill_value=fv)
             varmskID = ncid.createVariable('bdy_msk', 'f4', ('y', 'x', ), fill_value=fv)
         else:
diff --git a/pynemo/nemo_bdy_ncpop.py b/pynemo/nemo_bdy_ncpop.py
index 1f55c5dd98fa61ad81ce692a5327eedd06b7cd27..7e67dc2f491e361bbf32cbe75b42e2becf544ea8 100644
--- a/pynemo/nemo_bdy_ncpop.py
+++ b/pynemo/nemo_bdy_ncpop.py
@@ -18,7 +18,7 @@ def write_data_to_file(filename, variable_name, data):
     ncid = Dataset(filename, 'a', clobber=False, format='NETCDF4')
     count = data.shape
 
-    three_dim_variables = ['votemper', 'vosaline', 'N1p', 'N3n', 'N5s','thetao','so']
+    three_dim_variables = ['votemper', 'vosaline', 'N1p', 'N3n', 'N5s','thetao','so','uo','vo']
     two_dim_variables = ['sossheig', 'vobtcrtx', 'vobtcrty', 'iicethic', 'ileadfra', 'isnowthi','zos']
 
     if variable_name in three_dim_variables:
diff --git a/pynemo/profile.py b/pynemo/profile.py
index 91108281007c53e95ea68f683c2414d12aeca733..ab394b1a77632da39171b895c342b8a287e33268 100644
--- a/pynemo/profile.py
+++ b/pynemo/profile.py
@@ -490,8 +490,8 @@ def process_bdy(setup_filepath=0, mask_gui=False):
     
     emap = {}
     grd  = [  't',  'u',  'v']
-    pair = [ None, 'uv', 'uv'] # TODO: devolve this to the namelist?
-    
+    #pair = [ None, 'uv', 'uv'] # TODO: devolve this to the namelist?
+    pair = [None, None, None]
     # TODO: The following is a temporary stop gap to assign variables for both CMEMS downloads
     #  and existing variable names. In future we need a slicker way of determining the variables to extract.
     # Perhaps by scraping the .ncml file - this way biogeochemical tracers
diff --git a/unit_tests/namelist_unit_test_offset.bdy b/unit_tests/namelist_unit_test_offset.bdy
index b11b850c0cf00705fd82741e8704412ad45f4607..10d164b27b40d293c12864107655ab9091d3a728 100644
--- a/unit_tests/namelist_unit_test_offset.bdy
+++ b/unit_tests/namelist_unit_test_offset.bdy
@@ -30,17 +30,17 @@
 !------------------------------------------------------------------------------
 !  grid information 
 !------------------------------------------------------------------------------
-   sn_src_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_src_hgr_zps.nc'
-   sn_src_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_src_zgr_zps.nc'
-   sn_dst_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_offset_dst_hgr_zps.nc'
-   sn_dst_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_offset_dst_zgr_zps.nc'
-   sn_src_msk = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/mask.nc'
-   sn_bathy   = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_offset_dst_bathy.nc'
+   sn_src_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_src_hgr_zps.nc'
+   sn_src_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_src_zgr_zps.nc'
+   sn_dst_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_offset_dst_hgr_zps.nc'
+   sn_dst_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_offset_dst_zgr_zps.nc'
+   sn_src_msk = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/mask.nc'
+   sn_bathy   = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_offset_dst_bathy.nc'
 
 !------------------------------------------------------------------------------
 !  I/O 
 !------------------------------------------------------------------------------
-   sn_src_dir = '/Users/thopri/Projects/PyNEMO/unit_tests/src_data_unit_tests.ncml' ! src_files/'
+   sn_src_dir = '/Users/thopri/Projects/PyNEMO/unit_tests/unit_tests.ncml' ! src_files/'
    sn_dst_dir = '/Users/thopri/Projects/PyNEMO/unit_tests/test_outputs'
    sn_fn      = 'unit_test_offset'             ! prefix for output files
    nn_fv      = -1e20                 !  set fill value for output files
@@ -61,7 +61,7 @@
     ln_mask_file   = .false.              !  =T : read mask from file
     cn_mask_file   = 'mask.nc'            !  name of mask file 
                                           !  (if ln_mask_file=.TRUE.)
-    ln_dyn2d       = .false.              !  boundary conditions for 
+    ln_dyn2d       = .true.              !  boundary conditions for
                                           !  barotropic fields
     ln_dyn3d       = .false.              !  boundary conditions for 
                                           !  baroclinic velocities
diff --git a/unit_tests/namelist_unit_test_orth.bdy b/unit_tests/namelist_unit_test_orth.bdy
index 35e9d9d55341b97355da1d81c0a6bc3db3a769d7..854265223c1d86151fc7aa2807b754e894398417 100644
--- a/unit_tests/namelist_unit_test_orth.bdy
+++ b/unit_tests/namelist_unit_test_orth.bdy
@@ -30,17 +30,17 @@
 !------------------------------------------------------------------------------
 !  grid information 
 !------------------------------------------------------------------------------
-   sn_src_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_src_hgr_zps.nc'
-   sn_src_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_src_zgr_zps.nc'
-   sn_dst_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_dst_hgr_zps.nc'
-   sn_dst_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_dst_zgr_zps.nc'
-   sn_src_msk = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/mask.nc'
-   sn_bathy   = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_dst_bathy.nc'
+   sn_src_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_src_hgr_zps.nc'
+   sn_src_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_src_zgr_zps.nc'
+   sn_dst_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_dst_hgr_zps.nc'
+   sn_dst_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_dst_zgr_zps.nc'
+   sn_src_msk = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/mask.nc'
+   sn_bathy   = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_dst_bathy.nc'
 
 !------------------------------------------------------------------------------
 !  I/O 
 !------------------------------------------------------------------------------
-   sn_src_dir = '/Users/thopri/Projects/PyNEMO/unit_tests/src_data_unit_tests.ncml' ! src_files/'
+   sn_src_dir = '/Users/thopri/Projects/PyNEMO/unit_tests/unit_tests.ncml' ! src_files/'
    sn_dst_dir = '/Users/thopri/Projects/PyNEMO/unit_tests/test_outputs'
    sn_fn      = 'unit_test_orth'             ! prefix for output files
    nn_fv      = -1e20                 !  set fill value for output files
@@ -61,9 +61,9 @@
     ln_mask_file   = .false.              !  =T : read mask from file
     cn_mask_file   = 'mask.nc'            !  name of mask file 
                                           !  (if ln_mask_file=.TRUE.)
-    ln_dyn2d       = .false.              !  boundary conditions for
+    ln_dyn2d       = .true.              !  boundary conditions for
                                           !  barotropic fields
-    ln_dyn3d       = .false.              !  boundary conditions for 
+    ln_dyn3d       = .false.              !  boundary conditions for
                                           !  baroclinic velocities
     ln_tra         = .true.               !  boundary conditions for T and S
     ln_ice         = .false.              !  ice boundary condition   
diff --git a/unit_tests/namelist_unit_test_rotated.bdy b/unit_tests/namelist_unit_test_rotated.bdy
index 063affb9b3cc076275e46dfe34f0959c16ff0da5..f970101f3e7466628e7b8f35c179e0e55499592c 100644
--- a/unit_tests/namelist_unit_test_rotated.bdy
+++ b/unit_tests/namelist_unit_test_rotated.bdy
@@ -30,17 +30,17 @@
 !------------------------------------------------------------------------------
 !  grid information 
 !------------------------------------------------------------------------------
-   sn_src_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_src_hgr_zps.nc'
-   sn_src_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_src_zgr_zps.nc'
-   sn_dst_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_rot_dst_hgr_zps.nc'
-   sn_dst_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_rot_dst_zgr_zps.nc'
-   sn_src_msk = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/mask.nc'
-   sn_bathy   = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_rot_dst_bathy.nc'
+   sn_src_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_src_hgr_zps.nc'
+   sn_src_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_src_zgr_zps.nc'
+   sn_dst_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_rot_dst_hgr_zps.nc'
+   sn_dst_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_rot_dst_zgr_zps.nc'
+   sn_src_msk = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/mask.nc'
+   sn_bathy   = '/Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/test_rot_dst_bathy.nc'
 
 !------------------------------------------------------------------------------
 !  I/O 
 !------------------------------------------------------------------------------
-   sn_src_dir = '/Users/thopri/Projects/PyNEMO/unit_tests/src_data_unit_tests.ncml' ! src_files/'
+   sn_src_dir = '/Users/thopri/Projects/PyNEMO/unit_tests/unit_tests.ncml' ! src_files/'
    sn_dst_dir = '/Users/thopri/Projects/PyNEMO/unit_tests/test_outputs'
    sn_fn      = 'unit_test_rotated'             ! prefix for output files
    nn_fv      = -1e20                 !  set fill value for output files
@@ -61,7 +61,7 @@
     ln_mask_file   = .false.              !  =T : read mask from file
     cn_mask_file   = 'mask.nc'            !  name of mask file 
                                           !  (if ln_mask_file=.TRUE.)
-    ln_dyn2d       = .false.              !  boundary conditions for 
+    ln_dyn2d       = .true.              !  boundary conditions for
                                           !  barotropic fields
     ln_dyn3d       = .false.              !  boundary conditions for 
                                           !  baroclinic velocities
diff --git a/unit_tests/src_data_unit_tests.ncml b/unit_tests/src_data_unit_tests.ncml
deleted file mode 100644
index 3016c52e5a5d5da17889f49ac04430d7d57c5503..0000000000000000000000000000000000000000
--- a/unit_tests/src_data_unit_tests.ncml
+++ /dev/null
@@ -1,10 +0,0 @@
-<ns0:netcdf xmlns:ns0="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2" title="NEMO aggregation">
-  <ns0:aggregation type="union">
-    <ns0:netcdf>
-      <ns0:aggregation dimName="time" name="temperature" type="joinExisting">
-        <ns0:scan location="file://Users/thopri/Projects/PyNEMO/unit_tests/test_data/" regExp=".*T\.nc$" />
-      </ns0:aggregation>
-    </ns0:netcdf>
-  </ns0:aggregation>
-  <ns0:variable name="thetao" orgName="thetao" />
-</ns0:netcdf>
diff --git a/unit_tests/test_gen.py b/unit_tests/test_gen.py
index 0d6f98213109ab4f19f0650980e1fba1780b3046..9c27327be8ae3d3232b013f80c565d2c69fd35d9 100644
--- a/unit_tests/test_gen.py
+++ b/unit_tests/test_gen.py
@@ -21,8 +21,8 @@ def _main():
     max_dep = 100
     min_dep = 10
     z_end_dim = 1
-    h_fname = 'unit_tests/test_data/test_src_hgr_zps.nc'
-    z_fname = 'unit_tests/test_data/test_src_zgr_zps.nc'
+    h_fname = 'unit_tests/test_inputs/test_src_hgr_zps.nc'
+    z_fname = 'unit_tests/test_inputs/test_src_zgr_zps.nc'
     grid_h1 = gt.set_hgrid(dx,dy,jpi,jpj)
     grid_z1 = gt.set_zgrid(grid_h1,jpk,max_dep,min_dep,z_end_dim)
     write_coord_H = gt.write_coord_H(h_fname,grid_h1)
@@ -42,14 +42,14 @@ def _main():
     min_dep = 10
     z_end_dim = 1
     sf = 10
-    h_fname = 'unit_tests/test_data/test_dst_hgr_zps.nc'
-    z_fname = 'unit_tests/test_data/test_dst_zgr_zps.nc'
+    h_fname = 'unit_tests/test_inputs/test_dst_hgr_zps.nc'
+    z_fname = 'unit_tests/test_inputs/test_dst_zgr_zps.nc'
     grid_h2 = gt.set_hgrid(dx,dy,jpi,jpj,zoffx,zoffy,sf)
     grid_z2 = gt.set_zgrid(grid_h2,jpk,max_dep,min_dep,z_end_dim)
     write_coord_H = gt.write_coord_H(h_fname,grid_h2)
     write_coord_Z = gt.write_coord_Z(z_fname,grid_h2,grid_z2)
     # write bathy files (constant bathy)
-    bathy_fname = 'unit_tests/test_data/test_dst_bathy.nc'
+    bathy_fname = 'unit_tests/test_inputs/test_dst_bathy.nc'
     bathy = gt.write_bathy(bathy_fname,grid_h2,grid_z2)
     if write_coord_H + write_coord_Z + bathy == 0:
         print("Org child grid generation successful!")
@@ -60,8 +60,8 @@ def _main():
     origin = (8,8)
 
     # rotate grid
-    rot_h_fname = 'unit_tests/test_data/test_rot_dst_hgr_zps.nc'
-    rot_z_fname = 'unit_tests/test_data/test_rot_dst_zgr_zps.nc'
+    rot_h_fname = 'unit_tests/test_inputs/test_rot_dst_hgr_zps.nc'
+    rot_z_fname = 'unit_tests/test_inputs/test_rot_dst_zgr_zps.nc'
     grid_rot = grid_h2.copy()
     grid_rot['latt'], grid_rot['lont'] = gt.rotate_around_point(grid_h2['latt'],grid_h2['lont'],theta,origin)
     grid_rot['latu'], grid_rot['lonu'] = gt.rotate_around_point(grid_h2['latu'], grid_h2['lonu'], theta, origin)
@@ -70,7 +70,7 @@ def _main():
     write_coord_H = gt.write_coord_H(rot_h_fname,grid_rot)
     write_coord_Z = gt.write_coord_Z(rot_z_fname,grid_rot,grid_z2)
     # write bathy files (constant bathy)
-    bathy_fname = 'unit_tests/test_data/test_rot_dst_bathy.nc'
+    bathy_fname = 'unit_tests/test_inputs/test_rot_dst_bathy.nc'
     bathy = gt.write_bathy(bathy_fname,grid_rot,grid_z2)
     if write_coord_H + write_coord_Z + bathy == 0:
         print("Rotated child grid generation Successful!")
@@ -87,14 +87,14 @@ def _main():
     min_dep = 10
     z_end_dim = 1
     sf = 10
-    h_fname = 'unit_tests/test_data/test_offset_dst_hgr_zps.nc'
-    z_fname = 'unit_tests/test_data/test_offset_dst_zgr_zps.nc'
+    h_fname = 'unit_tests/test_inputs/test_offset_dst_hgr_zps.nc'
+    z_fname = 'unit_tests/test_inputs/test_offset_dst_zgr_zps.nc'
     grid_h3 = gt.set_hgrid(dx,dy,jpi,jpj,zoffx,zoffy,sf)
     grid_z3 = gt.set_zgrid(grid_h2,jpk,max_dep,min_dep,z_end_dim)
     write_coord_H = gt.write_coord_H(h_fname,grid_h3)
     write_coord_Z = gt.write_coord_Z(z_fname,grid_h3,grid_z3)
     # write bathy files (constant bathy)
-    bathy_fname = 'unit_tests/test_data/test_offset_dst_bathy.nc'
+    bathy_fname = 'unit_tests/test_inputs/test_offset_dst_bathy.nc'
     bathy = gt.write_bathy(bathy_fname,grid_h3,grid_z3)
     if write_coord_H + write_coord_Z + bathy == 0:
         print("Offset child grid gneration successful!")
@@ -104,7 +104,7 @@ def _main():
     #           grid_h3['lont'],grid_h1['latt'],grid_h1['lont'])
 
     # write boundary files (constant parameters)
-    out_fname = 'unit_tests/test_data/output_boundary' #drop file extension
+    out_fname = 'unit_tests/test_inputs/output_boundary' #drop file extension
     params_t = {'param1': {'name':'thetao','const_value':15.0,'longname':'temperature','units':'degreesC'},
               'param2': {'name':'so','const_value':35.0,'longname':'salinity','units':'PSU'},
               'param3': {'name': 'zos', 'const_value': 1.0, 'longname': 'sea surface height', 'units': 'metres'}
@@ -121,7 +121,7 @@ def _main():
         print('Boundary file generation successful!')
 
     #write_mask
-    mask_fname = 'unit_tests/test_data/mask.nc'
+    mask_fname = 'unit_tests/test_inputs/mask.nc'
     mask = gt.write_mask(mask_fname,grid_h1,grid_z1)
     if mask == 0:
         print('Mask file generation successful!')
diff --git a/unit_tests/unit_test.py b/unit_tests/unit_test.py
index 2f5a5c8d7ee91498fe3f1524576366646f6a4fa6..2f041ade946dbf909bfc3d28e6089d08534e599d 100644
--- a/unit_tests/unit_test.py
+++ b/unit_tests/unit_test.py
@@ -34,7 +34,9 @@ for n in namelist_files:
 
 # perform tests
 def test_temp():
-    test_files = glob.glob('unit_tests/test_outputs/unit_test*')
+    test_files = glob.glob('unit_tests/test_outputs/*bdyT*')
+    if len(test_files) == 0:
+        raise Exception('DONT PANIC: no temperature test files found')
     for t in test_files:
         results = Dataset(t) # open results
         temp = results['thetao'][:]
@@ -45,7 +47,9 @@ def test_temp():
         assert abs(temp_[temp_ != 0.0].min() - 15) <= 0.001
 
 def test_salinty():
-    test_files = glob.glob('unit_tests/test_outputs/unit_test*')
+    test_files = glob.glob('unit_tests/test_outputs/*bdyT*')
+    if len(test_files) == 0:
+        raise Exception('DONT PANIC: no salinity test files found')
     for t in test_files:
         results = Dataset(t)  # open results
         sal = results['so'][:]
@@ -55,18 +59,57 @@ def test_salinty():
         assert abs(sal_[sal_ != 0.0].max() - 35) <= 0.001
         assert abs(sal_[sal_ != 0.0].min() - 35) <= 0.001
 
+def test_ssh():
+    test_files = glob.glob('unit_tests/test_outputs/*bdyT*')
+    if len(test_files) == 0:
+        raise Exception('DONT PANIC: no SSH test files found')
+    for t in test_files:
+        results = Dataset(t)  # open results
+        ssh = results['zos'][:]
+        results.close()
+        ssh_ = np.ma.masked_array(ssh,ssh == -32767.0)
+        assert abs(ssh_[ssh_!=0.0].mean() - 1.0) <= 0.001
+        assert abs(ssh_[ssh_ != 0.0].max() - 1.0) <= 0.001
+        assert abs(ssh_[ssh_ != 0.0].min() - 1.0) <= 0.001
+
+def test_U():
+    test_files = glob.glob('unit_tests/test_outputs/*bdyU*')
+    if len(test_files) == 0:
+        raise Exception('DONT PANIC: no U current test files found')
+    for t in test_files:
+        results = Dataset(t)  # open results
+        U = results['uo'][:]
+        results.close()
+        U_ = np.ma.masked_array(U,U == -32767.0)
+        assert abs(U_[U_!=0.0].mean() - 0.5) <= 0.001
+        assert abs(U_[U_ != 0.0].max() - 0.5) <= 0.001
+        assert abs(U_[U_ != 0.0].min() - 0.5) <= 0.001
+
+def test_V():
+    test_files = glob.glob('unit_tests/test_outputs/*bdyV*')
+    if len(test_files) == 0:
+        raise Exception('DONT PANIC: no V current test files found')
+    for t in test_files:
+        results = Dataset(t)  # open results
+        V = results['vo'][:]
+        results.close()
+        V_ = np.ma.masked_array(V,V == -32767.0)
+        assert abs(V_[V_!=0.0].mean() - 0.5) <= 0.001
+        assert abs(V_[V_ != 0.0].max() - 0.5) <= 0.001
+        assert abs(V_[V_ != 0.0].min() - 0.5) <= 0.001
+
 # clean up test I/O
-files = glob.glob('unit_tests/test_outputs/*')
-for f in files:
-    os.remove(f)
-files = glob.glob('unit_tests/test_outputs/*')
-if len(files) != 0:
-    raise Exception('DONT PANIC: output folder not cleaned')
+def test_rm_out():
+    files = glob.glob('unit_tests/test_outputs/*')
+    for f in files:
+        os.remove(f)
+    files = glob.glob('unit_tests/test_outputs/*')
+    assert len(files) == 0
 
-files = glob.glob('unit_tests/test_data/*')
-for f in files:
-    os.remove(f)
-files = glob.glob('unit_tests/test_data/*')
-if len(files) != 0:
-    raise Exception('DONT PANIC: input folder not cleaned')
 
+def test_rm_in():
+    files = glob.glob('unit_tests/test_inputs/*')
+    for f in files:
+        os.remove(f)
+    files = glob.glob('unit_tests/test_inputs/*')
+    assert len(files) == 0
diff --git a/unit_tests/unit_tests.ncml b/unit_tests/unit_tests.ncml
new file mode 100644
index 0000000000000000000000000000000000000000..d376a10e7d9f6a9a1ae4511f74b400c994515b50
--- /dev/null
+++ b/unit_tests/unit_tests.ncml
@@ -0,0 +1,28 @@
+<ns0:netcdf xmlns:ns0="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2" title="NEMO aggregation">
+  <ns0:aggregation type="union">
+    <ns0:netcdf>
+      <ns0:aggregation dimName="time" name="temperature" type="joinExisting">
+        <ns0:scan location="file://Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/" regExp=".*T\.nc$" />
+      </ns0:aggregation>
+    </ns0:netcdf>
+    <ns0:netcdf>
+      <ns0:aggregation dimName="time" name="zonal_velocity" type="joinExisting">
+        <ns0:scan location="file://Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/" regExp=".*U\.nc$" />
+      </ns0:aggregation>
+    </ns0:netcdf>
+    <ns0:netcdf>
+      <ns0:aggregation dimName="time" name="meridian_velocity" type="joinExisting">
+        <ns0:scan location="file://Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/" regExp=".*V\.nc$" />
+      </ns0:aggregation>
+    </ns0:netcdf>
+    <ns0:netcdf>
+      <ns0:aggregation dimName="time" name="sea_surface_height" type="joinExisting">
+        <ns0:scan location="file://Users/thopri/Projects/PyNEMO/unit_tests/test_inputs/" regExp=".*T\.nc$" />
+      </ns0:aggregation>
+    </ns0:netcdf>
+  </ns0:aggregation>
+  <ns0:variable name="thetao" orgName="thetao" />
+  <ns0:variable name="uo" orgName="uo" />
+  <ns0:variable name="vo" orgName="vo" />
+  <ns0:variable name="zos" orgName="zos" />
+</ns0:netcdf>