diff --git a/Dockerfile b/Dockerfile
index b2dadfeb6039f839d21aa9f3a82744965728ea29..795026ffe2874473f71de2c0415292b23e403c7c 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -185,7 +185,6 @@ FROM ubuntu:20.04 as runner
 RUN mkdir /src
 COPY --from=builder /usr/local /usr/local
 RUN chmod +rx /usr/local/bin/*
-COPY --from=builder /src/wsclean /opt/
 
 SHELL ["/bin/bash", "-c"]
 
diff --git a/cluster.py b/cluster.py
index 898171aabfe773c5085d172a9753bad6a002e9c3..9b47100ec1b639d7d6ded6b3ba502c03b84451a3 100755
--- a/cluster.py
+++ b/cluster.py
@@ -605,21 +605,12 @@ def voronoi_clustering(fig, ax, df, wcs, resid_data, nbright, nclusters,
 
     # logging.info('Checking {} brightest model components'.format(nbright))
 
-    bright_df = df.sort_values('I')[::-1][:nbright][['ra', 'dec', 'I']]
-    # csnrs = []
-    # cfluxes = []
-    # cmeasures = [] # the main value to clusterize by
-    # cellipse_params = []
-
-    # fmin, fmax = min(a.I), max(a.I)
-    rms = mad(resid_data)
-    # resid_mean = np.mean(resid_data)
-
+    bright_df = df.sort_values('I')[::-1][['ra', 'dec', 'I']]
+    # rms = mad(resid_data)
     logging.info('Getting measures for the potential clusters...')
     clusters = []
     clusters_centers = [] #
 
-
     for ra, dec, flux in bright_df.values:
         c = SkyCoord(ra, dec, unit='deg')
         px, py = np.round(wcs.all_world2pix(c.ra, c.dec, 0)).astype(int)
@@ -627,20 +618,18 @@ def voronoi_clustering(fig, ax, df, wcs, resid_data, nbright, nclusters,
         if (abs(px-resid_data.shape[1]) < boxsize) or (abs(py-resid_data.shape[0]) < boxsize):
             # logging.debug('Skipping the edge source')
             continue
-# Check if the component already in a cluster
+# Check if the component is nearby
         if clusters and any([c.separation(_).arcmin<same_source_radius for _ in clusters]):
             continue
 
         small_resid = resid_data[py-boxsize:py+boxsize, px-boxsize:px+boxsize]
         ellipse_mean, ecc, amaj, numpix = ellipses_coh(small_resid, amin=20, amax=boxsize-1, dr=1.0)
 
-        if abs(ellipse_mean/rms) > 1.4:
-             # rect = plt.Rectangle((px-boxsize, py-boxsize), 2*boxsize, 2*boxsize,
-             #                      lw=1, color='k', fc='none', alpha=0.3)
-             # ax.add_artist(rect)
-             clusters_centers.append([ra, dec])
-             clusters.append(c)
-             print(ra, dec)
+# Uncomment the condition to search for artifacts:
+        # if abs(ellipse_mean/rms) > 1.4:
+        clusters_centers.append([ra, dec])
+        clusters.append(c)
+        print(ra, dec)
 
         if (isinstance(nclusters, int)) and (len(clusters_centers) >= nclusters):
             logging.debug('Max cluster number reached. Breaking...')
@@ -650,12 +639,6 @@ def voronoi_clustering(fig, ax, df, wcs, resid_data, nbright, nclusters,
         logging.warning('Decreasing number of clusters')
         nclusters = len(clusters_centers)
 
-    # X = np.vstack([df.ra, df.dec]).T
-    # kmeans = KMeans(n_clusters=nclusters, init=clusters_centers, n_init=1, random_state=0, max_iter=1)
-    # kmeans.fit(X)
-    # y_kmeans = kmeans.predict(X) # cluster index for each observation
-    # ax.scatter(clusters_centers[:, 0], clusters_centers[:, 1], c='black', s=100, alpha=0.5, transform=ax.get_transform("world"))
-
     vor = Voronoi(np.array(clusters_centers))
     voronoi_plot_2d_world(vor, ax=ax, show_vertices=False)
 
@@ -865,22 +848,22 @@ def write_ds9(fname, h5, image, points=None):
     xmax = imheader['NAXIS1']
     ymin = 0
     ymax = imheader['NAXIS2']
- 
+
     # To cut the Voronoi tessellation on the bounding box, we need
     # a "circumscribing circle"
     dist_pix = np.sqrt((xmax - xmin) ** 2 + (ymax - ymin) ** 2)
-    
-    
+
+
     # load in the directions from the H5
     sourcedir = read_dir_fromh5(h5)
-    
+
     # make ra and dec arrays and coordinates c
     ralist = sourcedir[:, 0]
     declist = sourcedir[:, 1]
     c = SkyCoord(ra=ralist * u.rad, dec=declist * u.rad)
-    
+
     # convert from ra,dec to x,y pixel
-    x, y = w.wcs_world2pix(c.ra.degree, c.dec.degree, 1)    
+    x, y = w.wcs_world2pix(c.ra.degree, c.dec.degree, 1)
 
     bbox = Polygon([(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax)])
     polygons, points = tessellate(x, y, w, dist_pix, bbox, plot_tessellation=False)
diff --git a/imcal.py b/imcal.py
index 2ceb05c7c45d10a8882321f17899696ba2bd4432..f929ac8e0d1e73369eefe37b039fc7d54c56733f 100755
--- a/imcal.py
+++ b/imcal.py
@@ -40,7 +40,6 @@ _MAX_TIME = 1 * 3600 # SECONDS
 _MAX_POOL = _MAX_TIME // _POOL_TIME
 
 def execute_binary(binary, args, ):
-    
     command = [f'{binary}'] + args
     logging.debug('executing %s', ','.join(command))
     dppp_process = subprocess.Popen(command)
@@ -51,7 +50,7 @@ def execute_binary(binary, args, ):
             return return_code
         except TimeoutExpired as e:
             logging.debug('DPPP process %s still running', dppp_process.pid)
-            continue 
+            continue
 
 
 def fft_psf(bmaj, bmin, bpa, size=3073):
@@ -65,6 +64,7 @@ def fft_psf(bmaj, bmin, bpa, size=3073):
     fkern = fkern.array
     return fkern
 
+
 def reconvolve_gaussian_kernel(img, old_maj, old_min, old_pa, new_maj, new_min, new_pa):
     """
     convolve image with a gaussian kernel without FFTing it
@@ -123,7 +123,7 @@ def modify_filename(fname, string, ext=None):
     return fbase + string + fext
 
 
-def wsclean(msin, wsclean_bin='wsclean', datacolumn='DATA', outname=None, pixelsize=3, imagesize=3072, mgain=0.8, 
+def wsclean(msin, wsclean_bin='wsclean', datacolumn='DATA', outname=None, pixelsize=3, imagesize=3072, mgain=0.8,
             multifreq=0, autothresh=0.3,
             automask=3, niter=1000000, multiscale=False, save_source_list=True,
             clearfiles=True, clip_model_level=None,
@@ -178,12 +178,13 @@ def smoothImage(imgfits, psf=30, out=None) :
     return fits_reconvolve_psf(imgfits, Beam(psf*u.arcsec), out=out)
 
 
-def create_mask(imgfits, residfits, clipval, outname='mask.fits', ):
+def create_mask(imgfits, residfits, clipval, outname='mask.fits'):
     """
     Create mask using Tom's code (e-mail on 1 Jul 2021)
     """
-    cmd1 = f'makeNoiseMapFitsLow {imgfits} {residfits} noise.fits noiseMap.fits'
-    cmd2 = f'makeMaskFits noiseMap.fits {outname} {clipval}'
+    outbase = os.path.splitext(imgfits)[0]
+    cmd1 = f'makeNoiseMapFitsLow {imgfits} {residfits} {outbase}_noise.fits {outbase}_noiseMap.fits'
+    cmd2 = f'makeMaskFits {outbase}_noiseMap.fits {outname} {clipval}'
     logging.debug("Running command: %s", cmd1)
     subprocess.call(cmd1, shell=True)
     logging.debug("Running command: %s", cmd2)
@@ -191,42 +192,50 @@ def create_mask(imgfits, residfits, clipval, outname='mask.fits', ):
     return outname
 
 
-def makeNoiseImage(imgfits, residfits, low=False, ) :
+def makeNoiseImage(imgfits, residfits, low=False) :
     """
     Create mask using Tom's code (e-mail on 1 Jul 2021)
     """
+    # if outbase is None:
+    outbase = os.path.splitext(imgfits)[0]
     if low:
-      cmd = f'makeNoiseMapFitsLow {imgfits} {residfits} noiseLow.fits noiseMapLow.fits'
+        img1, img2 = f'{outbase}_noiseLow.fits', f'{outbase}_noiseMapLow.fits'
+        cmd = f'makeNoiseMapFitsLow {imgfits} {residfits} {img1} {img2}'
     else :
-      cmd = f'makeNoiseMapFits {imgfits} {residfits} noise.fits noiseMap.fits'
-      
+        img1, img2 = f'{outbase}_noise.fits', f'{outbase}_noiseMap.fits'
+        cmd = f'makeNoiseMapFits {imgfits} {residfits} {img1} {img2}'
     logging.debug("Running command: %s", cmd)
     subprocess.call(cmd, shell=True)
-    return
+    return img2
 
 
-def makeCombMask(ima1='noiseMap.fits', ima2='noiseMapLow.fits',
-		 clip1=5, clip2=7, outname='mask.fits', ) :
+def makeCombMask(img1, img2, clip1=5, clip2=7, outname=None) :
     """
     Create mask using Tom's code (e-mail on 1 Jul 2021)
     """
-    cmd = f'makeCombMaskFits {ima1} {ima2} {outname} {clip1} {clip2}'
-    
+    if outname is None:
+        outname = os.path.splitext(img1)[0] + '_mask.fits'
+    cmd = f'makeCombMaskFits {img1} {img2} {outname} {clip1} {clip2}'
     logging.debug("Running command: %s", cmd)
     subprocess.call(cmd, shell=True)
-    return
+    return outname
 
 
 def get_image_ra_dec_min_max(msin):
     """
     Determine image center coords, min and max values for msin
     """
-    cmd = f'wsclean -niter 0 -size 3072 3072 -scale 3arcsec -use-wgridder {msin}'
-    logging.debug('Running command: %s', cmd)
-    subprocess.call(cmd, shell=True)
-    data = fits.getdata('wsclean-image.fits')
-    header = fits.getheader('wsclean-image.fits')
-    return header['CRVAL1'], header['CRVAL2'], np.nanmin(data), np.nanmax(data)
+    outbase = os.path.splitext(msin.rstrip('/'))[0]+'-wsclean'
+    outname = outbase+'-image.fits'
+    cmd = f'wsclean -name {outbase} -niter 0 -size 3072 3072 -scale 3arcsec -use-wgridder {msin}'
+    if os.path.exists(outname):
+        logging.debug('Image exists. Skipping cleaning...')
+    else:
+        logging.debug('Running command: %s', cmd)
+        subprocess.call(cmd, shell=True)
+    data = fits.getdata(outname)
+    header = fits.getheader(outname)
+    return outname, header['CRVAL1'], header['CRVAL2'], np.nanmin(data), np.nanmax(data)
 
 
 def makesourcedb(modelfile, out=None, ):
@@ -255,9 +264,7 @@ def render(bkgr, model, out=None, ):
     return out
 
 
-        
 def execute_dppp(args, ):
-    
     command = ['DP3'] + args
     logging.debug('executing %s', ','.join(command))
     dppp_process = subprocess.Popen(command)
@@ -298,6 +305,25 @@ def split_ms(msin_path, startchan, nchan=0, msout_path='', ):
     return msout_path
 
 
+def preflag(msin, msout=None, **kwargs):
+    """
+    preflag data using DP3 preflag module
+    """
+    if not any(kwargs.values()):
+        logging.debug('No preflag options specified. Skipping...')
+        return msin
+    msout = msout or '.'
+    command_args = ['steps=[preflag]',
+                    f'msin={msin}',
+                    f'msout={msout}'] + ['preflag.'+'='.join(_) for _ in kwargs.items() if _[1] is not None]
+    logging.info('Flagging data (%s)', command_args)
+    return_code = execute_dppp(command_args)
+    logging.debug('Preflag of %s returned status code %s', msin, return_code)
+    check_return_code(return_code)
+    if msout == '.': msout = msin
+    return msout
+
+
 def dical(msin, srcdb, msout=None, h5out=None, solint=1, startchan=0, split_nchan=0,
           mode='phaseonly', cal_nchan=0, uvlambdamin=500, ):
     """ direction independent calibration with DPPP """
@@ -426,7 +452,7 @@ def view_sols(h5param, outname=None):
             fig2 = ax2 = None
             logging.error('No phase solutions found')
     # return fig1, ax1, fig2, ax2
-    
+
 
 def remove_model_components_below_level(model, level=0.0, out=None):
     """
@@ -465,9 +491,9 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
     with open(cfgfile) as f:
         cfg = yaml.safe_load(f)
 
-    
+
     if steps == 'all':
-        steps = ['nvss', 'mask', 'dical', 'ddcal']
+        steps = ['nvss', 'preflag', 'mask', 'dical', 'ddcal']
     else:
         steps = steps.split(',')
 
@@ -480,6 +506,7 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
 
     ms_split = msbase + '_splt.MS'
 
+
     img0 = outbase + '_0'
     img1 = outbase + '_1'
     img2 = outbase + '_2'
@@ -487,8 +514,8 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
     img_dical = outbase + '-dical'
     img_ddsub_1 = outbase + '-ddsub-1'
     img_ddsub_2 = outbase + '-ddsub-2'
-    img_ddcal_1 = outbase + '-ddcal'
-    img_ddcal_2 = outbase + '-ddcal'
+    img_ddcal_1 = outbase + '-ddcal-1'
+    img_ddcal_2 = outbase + '-ddcal-2'
 
     mask0 = outbase + '-mask0.fits'
     mask1 = outbase + '-mask1.fits'
@@ -514,16 +541,21 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
     h5_dd = outbase + '_ddcal.h5'
 
 
-    
+
     if not force and os.path.exists(img_ddcal_2+'-image.fits'):
         logging.info('The final image exists. Exiting...')
         return 0
 
 # get image parameters
-    img_ra, img_dec, img_min, img_max = get_image_ra_dec_min_max(msin)
-    
+    # if 'init' in steps:
+    # if not force and os.path.exists(outbase+'-wsclean-image.fits'):
+    initial_img, img_ra, img_dec, img_min, img_max = get_image_ra_dec_min_max(msin)
+    logging.info('Image: %s', initial_img)
+    logging.info('Image RA, DEC: %s, %s', img_ra, img_dec)
+    logging.info('Image Min, Max: %s, %s', img_min, img_max)
+
     if 'nvss' in steps and cfg['nvss']:
-       nvss_model = nvss_cutout('wsclean-image.fits', nvsscat='/opt/nvss.csv.zip', clip=cfg['nvsscal']['clip_model'])
+       nvss_model = nvss_cutout(initial_img, nvsscat='/opt/nvss.csv.zip', clip=cfg['nvsscal']['clip_model'])
        if (not os.path.exists(ms_split)) and (cfg['split']['startchan'] or cfg['split']['nchan']):
            ms_split = split_ms(msin, msout_path=ms_split, **cfg['split'])
 
@@ -532,8 +564,13 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
        dical(ms_split, nvssMod, msout=dical0, h5out=h5_0, **cfg['nvsscal'])
        view_sols(h5_0, outname=msbase+'_sols_dical0')
     else :
-       # if (not os.path.exists(ms_split)) and (cfg['split']['startchan'] or cfg['split']['nchan']):
-       ms_split = split_ms(msin, msout_path=dical0, **cfg['split'])
+       if (not os.path.exists(ms_split)) and (cfg['split']['startchan'] or cfg['split']['nchan']):
+           dical0 = split_ms(msin, msout_path=dical0, **cfg['split'])
+       else:
+           dical0 = msin
+
+    if 'preflag' in steps:
+        dical0 = preflag(ms_split, msout=dical0, **cfg['preflag'])
 
     if 'mask' in steps:
         if not force and (os.path.exists(img0 +'-image.fits') or (os.path.exists(img0 +'-MFS-image.fits'))):
@@ -544,7 +581,7 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
             wsclean(dical0, outname=img0, automask=None, save_source_list=False, multifreq=False, mgain=None,
                     kwstring=f'-threshold {threshold}')
             create_mask(img0 +'-image.fits', img0 +'-residual.fits', clipval=10, outname=mask0, )
-    
+
     if 'dical' in steps:
 # clean1
         if not force and (os.path.exists(img1 +'-image.fits') or (os.path.exists(img1 +'-MFS-image.fits'))):
@@ -552,7 +589,7 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         else:
             wsclean(dical0, fitsmask=mask0, outname=img1, **cfg['clean1']) # fast shallow clean
             makesourcedb(img1+'-sources.txt', out=model1)
-    
+
 # dical1
         if not force and os.path.exists(dical1):
             logging.debug('dical/dical1 step: MS exists, , use --f to overwrite...')
@@ -565,10 +602,10 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         else:
             wsclean(dical1, fitsmask=mask0, outname=img2, **cfg['clean2'])
             smoothImage(img2+'-residual.fits')
-            makeNoiseImage(img2 +'-image.fits', img2 +'-residual.fits', )
-            makeNoiseImage(img2 +'-residual-smooth.fits', img2 +'-residual.fits', low=True, )
-            makeCombMask(outname=mask1, clip1=7, clip2=15, )
-    
+            i1 = makeNoiseImage(img2 +'-image.fits', img2 +'-residual.fits', )
+            i2 = makeNoiseImage(img2 +'-residual-smooth.fits', img2 +'-residual.fits', low=True, )
+            makeCombMask(i1, i2, clip1=7, clip2=15, outname=mask1, )
+
             makesourcedb(img2+'-sources.txt', out=model2, )
 
 # dical2
@@ -583,10 +620,10 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         else:
             wsclean(dical2, fitsmask=mask1, outname=img3, **cfg['clean3'])
             smoothImage(img3+'-residual.fits')
-            makeNoiseImage(img3 +'-image.fits', img3 +'-residual.fits', )
-            makeNoiseImage(img3 +'-residual-smooth.fits', img3 +'-residual.fits', low=True, )
-            makeCombMask(outname=mask2, clip1=5, clip2=10)
-    
+            i1 = makeNoiseImage(img3 +'-image.fits', img3 +'-residual.fits', )
+            i2 = makeNoiseImage(img3 +'-residual-smooth.fits', img3 +'-residual.fits', low=True, )
+            makeCombMask(i1, i2, clip1=5, clip2=10, outname=mask2,)
+
             makesourcedb(img3+'-sources.txt', out=model3)
 
 # dical3
@@ -602,9 +639,9 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         else:
             wsclean(dical3, fitsmask=mask2, outname=img_dical,  **cfg['clean4'])
             smoothImage(img_dical+'-residual.fits')
-            makeNoiseImage(img_dical +'-image.fits', img_dical +'-residual.fits', )
-            makeNoiseImage(img_dical +'-residual-smooth.fits', img_dical +'-residual.fits',low=True, )
-            makeCombMask(outname=mask3, clip1=5, clip2=7, )
+            i1 = makeNoiseImage(img_dical +'-image.fits', img_dical +'-residual.fits', )
+            i2 = makeNoiseImage(img_dical +'-residual-smooth.fits', img_dical +'-residual.fits',low=True, )
+            makeCombMask(i1, i2, clip1=5, clip2=7, outname=mask3)
 
     if 'ddcal' in steps:
 # Cluster
@@ -627,35 +664,36 @@ def main(msin, steps='all', outbase=None, cfgfile='imcal.yml', force=False):
         if not force and os.path.exists(img_ddsub_1+'-image.fits'):
             pass
         else:
-            wsclean(ddsub, fitsmask=mask3,outname=img_ddsub_1, **cfg['clean5'])
+            wsclean(ddsub, fitsmask=mask3, outname=img_ddsub_1, **cfg['clean5'])
 #TAO        wsclean(ddsub,outname=img_ddsub, **cfg['clean5'])
 
         aomodel = bbs2model(img_dical+'-sources.txt', img_dical+'-model.ao', )
-    
+
         render(img_ddsub_1+'-image.fits', aomodel, out=img_ddcal_1+'-image.fits')
-    
-        smoothImage(img_ddsub_1+'-image.fits')
-        makeNoiseImage(img_ddcal_1 +'-image.fits', img_ddsub_1 +'-residual.fits', )
-        makeNoiseImage(img_ddcal_1 +'-image-smooth.fits', img_ddsub_1 +'-residual.fits',low=True, )
-        makeCombMask(outname=mask4, clip1=3.5, clip2=5, )
+
+        smoothImage(img_ddcal_1+'-image.fits')
+        i1 = makeNoiseImage(img_ddcal_1 +'-image.fits', img_ddsub_1 +'-residual.fits', )
+        i2 = makeNoiseImage(img_ddcal_1 +'-image-smooth.fits', img_ddsub_1 +'-residual.fits',low=True, )
+        makeCombMask(i1, i2, clip1=3.5, clip2=5, outname=mask4,)
 
         if not force and os.path.exists(img_ddsub_2+'-image.fits'):
             pass
         else:
             wsclean(ddsub, fitsmask=mask4, outname=img_ddsub_2, **cfg['clean5'])
 #TAO        wsclean(ddsub,outname=img_ddsub, **cfg['clean5'])
-    
+
         aomodel = bbs2model(img_dical+'-sources.txt', img_dical+'-model.ao', )
         render(img_ddsub_2+'-image.fits', aomodel, out=img_ddcal_2+'-image.fits', )
 
 # test facet imaging:
-    # ds9_file = 'ddfacets.reg'
-    # ddvis = outbase + '_ddvis.MS'
-    # h5_ddvis = 'ddsols.h5'
-    # clustered_sdb = img_dical+'-clustered.sourcedb'
-    # # ddvis = ddecal(dical3, clustered_sdb, msout=ddvis, subtract=False, h5out=h5_ddvis, **cfg['ddcal'])
-    # # write_ds9(ds9_file, h5_ddvis, img_ddcal+'-image.fits')
-    # wsclean(ddvis, fitsmask=mask3, save_source_list=False, outname='img-facet', **cfg['clean6'],)
+    if 'facet' in steps:
+        ds9_file = 'ddfacets.reg'
+        ddvis = outbase + '_ddvis.MS'
+        h5_ddvis = 'ddsols.h5'
+        clustered_sdb = img_dical+'-clustered.sourcedb'
+        # ddvis = ddecal(dical3, clustered_sdb, msout=ddvis, subtract=False, h5out=h5_ddvis, **cfg['ddcal'])
+        # write_ds9(ds9_file, h5_ddvis, img_ddcal_2+'-image.fits')
+        wsclean(ddvis, fitsmask=mask3, save_source_list=False, outname='img-facet', **cfg['facet_clean'],)
 
 
     return 0
@@ -673,13 +711,13 @@ if __name__ == "__main__":
                         dest='configfile', help='Config file', type=str)
     parser.add_argument('-o', '--outbase', default=None, help='output prefix', type=str)
     parser.add_argument('-s', '--steps', default='all', help='steps to run. Example: "nvss,mask,dical,ddcal"', type=str)
-    parser.add_argument('-f', '--force', action='store_false', help='Overwrite the existing files')
+    parser.add_argument('-f', '--force', action='store_true', help='Overwrite the existing files')
 
     args = parser.parse_args()
     configfile = args.configfile or \
         os.path.join(os.path.dirname(os.path.realpath(__file__)), 'imcal.yml')
     # msin = args.msin
-    main(args.msin, outbase=args.outbase, steps=args.steps, cfgfile=configfile)
+    main(args.msin, outbase=args.outbase, steps=args.steps, cfgfile=configfile, force=args.force)
     # extime = Time.now() - t0
     # print("Execution time: {:.1f} min".format(extime.to("minute").value))
 
diff --git a/imcal.yml b/imcal.yml
index 7ae17e911dbac5fa3e3d5019cdcc7ca6196431e5..bf0a4e8a4ac88e3b8ea60463b11244102f10e0cc 100644
--- a/imcal.yml
+++ b/imcal.yml
@@ -23,7 +23,14 @@ clean0: # initial clean
     max_over_thresh: 250 # the threshold for initial CLEAN is set to image_max/max_over_thresh
 # TODO add fixed threshold of 100 uJy
 
-############################## SPLIT ##########################################
+############################## PREFLAG ##########################################
+
+preflag: # See DP3 steps.preflag
+    baseline: # e.g. '[[RTB,*]]'
+    reltime: # e.g. '10:20+-20m'
+    chan: # null
+    blmin: # in meters
+
 clean1: # wsclean setup
     imagesize: 3072
     pixelsize: 3
@@ -86,7 +93,7 @@ clean4:
 
 ####################### CLUSTERING #######################
 cluster:
-    nbright: 80 # number of brightest clean components (CC) to check for artefacts
+    nbright: 100 # number of brightest clean components (CC) to check for artefacts
     boxsize: 250 # the boxsize around CC in pixels where to check for artefacts
     nclusters: 10 # number of clusters ('auto' -- to set automatically)
 # the following is only for 'auto' and 'manual' mathods:
@@ -114,16 +121,16 @@ clean5:
     clip_model_level: null
     kwstring: '-use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs -1.5'
 
-
-# clean6:
-#     imagesize: 3072
-#     pixelsize: 3
-#     multifreq: 6
-#     automask: 3.5
-#     autothresh: 0.5
-#     multiscale: True
-#     clip_model_level: null
-#     kwstring: '-facet-regions ddfacets.reg -apply-facet-solutions ddsols.h5 amplitude000,phase000 -use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs -1.5'
+### test for facet clean
+facet_clean:
+    imagesize: 3072
+    pixelsize: 3
+    multifreq: 6
+    automask: 3.5
+    autothresh: 0.5
+    multiscale: True
+    clip_model_level: null
+    kwstring: '-facet-regions ddfacets.reg -apply-facet-solutions ddsols.h5 amplitude000,phase000 -use-wgridder -parallel-deconvolution 1400 -parallel-gridding 8 -deconvolution-channels 3 -weight briggs -1.5'
 
 
 ### END