From a45eac37c0a3db1dfb2f14e7746719b569e8da74 Mon Sep 17 00:00:00 2001
From: alex <alex@tls-tautenburg.de>
Date: Mon, 31 Jul 2023 17:06:00 +0000
Subject: [PATCH] Add LINC LBA target imaging

---
 README.md                                  |   2 +-
 docs/source/running.rst                    |   2 +-
 docs/source/target.rst                     |   2 +
 scripts/blank_image_reg.py                 | 113 ++++
 scripts/download_skymodel_target.py        |   8 +-
 scripts/make_beam_reg.py                   |  95 +++
 scripts/make_mask.py                       |  70 +++
 scripts/sort_times_into_freqGroups.py      |   2 +-
 steps/applycal.cwl                         |   6 +
 steps/blank_image_reg.cwl                  |  70 +++
 steps/find_skymodel_target.cwl             |   9 +-
 steps/getFWHM.cwl                          |  71 +++
 steps/make_beam_reg.cwl                    |  55 ++
 steps/make_mask.cwl                        |  98 +++
 steps/sort_times_into_freqGroups.cwl       |  16 +
 steps/wsclean.cwl                          | 156 ++++-
 steps/wsclean_predict.cwl                  | 360 +++++++++++
 workflows/HBA_target.cwl                   |   8 +
 workflows/LBA_target.cwl                   |   9 +
 workflows/linc_target.cwl                  |  35 +-
 workflows/linc_target/corrupt_model.cwl    | 167 ++++++
 workflows/linc_target/finalize.cwl         |   2 +
 workflows/linc_target/gsmcal.cwl           |  76 +--
 workflows/linc_target/imaging_subtract.cwl | 666 +++++++++++++++++++++
 workflows/linc_target/prep.cwl             |  64 +-
 workflows/linc_target/selfcal_targ.cwl     |  63 ++
 26 files changed, 2167 insertions(+), 58 deletions(-)
 create mode 100755 scripts/blank_image_reg.py
 create mode 100755 scripts/make_beam_reg.py
 create mode 100755 scripts/make_mask.py
 create mode 100644 steps/blank_image_reg.cwl
 create mode 100644 steps/getFWHM.cwl
 create mode 100644 steps/make_beam_reg.cwl
 create mode 100644 steps/make_mask.cwl
 create mode 100644 steps/wsclean_predict.cwl
 create mode 100644 workflows/linc_target/corrupt_model.cwl
 create mode 100644 workflows/linc_target/imaging_subtract.cwl

diff --git a/README.md b/README.md
index f162eb5d..cf4f2f6d 100644
--- a/README.md
+++ b/README.md
@@ -32,7 +32,7 @@ The full documentation can be found at the [LINC webpage](https://linc.readthedo
 * [LofarStMan](https://github.com/lofar-astron/LofarStMan)
 * [Dysco](https://github.com/aroffringa/dysco.git) (v1.2 or later)
 * casacore
-* Python3 (including matplotlib, scipy, and astropy)
+* Python3 (including matplotlib, scipy, astropy)
 * [cwltool](https://github.com/common-workflow-language/cwltool) (3.1.20220406080846 or later) [toil-cwl-runner](https://github.com/DataBiosphere/toil) (5.7.0a1 or later)
 
 ### Installation
diff --git a/docs/source/running.rst b/docs/source/running.rst
index eecdc21f..dc5ccf8a 100644
--- a/docs/source/running.rst
+++ b/docs/source/running.rst
@@ -44,7 +44,7 @@ The following list provides the workflows to call in the command above for stand
 
 .. note::
     
-    The **LBA** target workflow is not (yet) available.
+    The **LBA** target workflow is still experimental and thus may not provide the expected results.
     
 If you have installed ``cwltool`` or ``toil`` locally on your system, **LINC** will pull automatically the right (u)Docker/Singularity image for you.
 
diff --git a/docs/source/target.rst b/docs/source/target.rst
index e2a9ef10..b04dfa6e 100644
--- a/docs/source/target.rst
+++ b/docs/source/target.rst
@@ -282,7 +282,9 @@ User-defined parameter configuration
 - ``updateweights``: update ``WEIGHT_SPECTRUM`` column in a way consistent with the weights being inverse proportional to the autocorrelations (default: ``true``)
 - ``use_target``: enable downloading of a target skymodel (default: ``true``)
 - ``skymodel_source``: choose the target skymodel from `TGSS ADR`_ or the new `Global Sky Model`_ (GSM) (default: ``TGSS``)
+- ``skymodel_fluxlimit``: limits the input skymodel to sources that exceed the given flux density limit in Jy (default: None for **HBA**, i.e. all sources of the catalogue will be kept, and 1.0 for **LBA**)
 - ``selfcal``: perform extensive self-calibration according to the `LiLF`_ scheme (recommended for **LBA** observations) (default: ``false``)
+- ``selfcal_region``: ds9-compatible region file to select the image regions used for the self-calibration
 
 A comprehensive explanation of the baseline selection syntax can be found `here`_.
 
diff --git a/scripts/blank_image_reg.py b/scripts/blank_image_reg.py
new file mode 100755
index 00000000..205de138
--- /dev/null
+++ b/scripts/blank_image_reg.py
@@ -0,0 +1,113 @@
+#!/usr/bin/env python
+
+# blank areas of a fits image given a region file
+
+# @author: Francesco de Gasperin, modified by Alexander Drabent
+
+def flatten(f, channel = 0, freqaxis = 0):
+    """
+    Flatten a fits file so that it becomes a 2D image. Return new header and data
+    """
+    from astropy import wcs
+    import numpy as np
+
+    naxis=f[0].header['NAXIS']
+    if (naxis < 2):
+        raise RadioError("Can\'t make map from this")
+    if (naxis == 2):
+        return f[0].header,f[0].data
+
+    w               = wcs.WCS(f[0].header)
+    wn              = wcs.WCS(naxis = 2)
+
+    wn.wcs.crpix[0] = w.wcs.crpix[0]
+    wn.wcs.crpix[1] = w.wcs.crpix[1]
+    wn.wcs.cdelt    = w.wcs.cdelt[0:2]
+    wn.wcs.crval    = w.wcs.crval[0:2]
+    wn.wcs.ctype[0] = w.wcs.ctype[0]
+    wn.wcs.ctype[1] = w.wcs.ctype[1]
+
+    header = wn.to_header()
+    header["NAXIS"] = 2
+    header["NAXIS1"] = f[0].header['NAXIS1']
+    header["NAXIS2"] = f[0].header['NAXIS2']
+    copy=('EQUINOX','EPOCH')
+    for k in copy:
+        r = f[0].header.get(k)
+        if (r):
+            header[k] = r
+
+    slicing = []
+    for i in range(naxis,0,-1):
+        if (i <= 2):
+            slicing.append(np.s_[:],)
+        elif (i == freqaxis):
+            slicing.append(channel)
+        else:
+            slicing.append(0)
+
+    # slice=(0,)*(naxis-2)+(np.s_[:],)*2
+    return header, f[0].data[tuple(slicing)]
+
+def main(input_image, region, outfile = None, inverse = False, blankval = 0., operation = "AND"):
+    """
+    Set to "blankval" all the pixels inside the given region
+    if inverse=True, set to "blankval" pixels outside region.
+    If a list of regions is provided the operation is applied to each region one after the other
+
+    input_image: fits file
+    region: ds9 region or list of regions
+    outfile: output name
+    inverse: reverse final *combined* mask
+    blankval: pixel value to set
+    operation: how to combine multiple regions with AND or OR
+    """
+
+    import astropy.io.fits as pyfits
+    import numpy as np
+    import pyregion
+
+
+    if outfile == None: outfile = input_image[0]
+    if not type(region) is list: region=[region]
+
+    # open fits
+    with pyfits.open(input_image[0]) as fits:
+        origshape    = fits[0].data.shape
+        header, data = flatten(fits)
+        sum_before   = np.sum(data)
+        if (operation == 'AND'):
+            total_mask = np.ones(shape = data.shape).astype(bool)
+        if (operation == 'OR'):
+            total_mask = np.zeros(shape = data.shape).astype(bool)
+        for this_region in region:
+            # extract mask
+            r    = pyregion.open(this_region)
+            mask = r.get_mask(header=header, shape=data.shape)
+            if (operation == 'AND'):
+                total_mask = total_mask & mask
+            if (operation == 'OR'):
+                total_mask = total_mask | mask
+        if (inverse):
+            total_mask = ~total_mask
+        data[total_mask] = blankval
+        # save fits
+        fits[0].data = data.reshape(origshape)
+        fits.writeto(outfile, overwrite=True)
+
+    print("%s: Blanking (%s): sum of values: %f -> %f" % (outfile, region, sum_before, np.sum(data)))
+
+if __name__=='__main__':
+    import argparse
+    parser = argparse.ArgumentParser(description='Set to "blankval" all the pixels inside the given region')
+    
+    parser.add_argument('input_image', type=str, nargs='+', help='Input FITS image')
+    parser.add_argument('region', type=str, nargs='+', help='Input region file')
+    parser.add_argument('--outfile', type=str, default=None, help='Name of output file, default: None')
+    parser.add_argument('--inverse', type=bool, default=False, help='Reverse final combined mask')
+    parser.add_argument('--blankval', type=float, default=0., help='Pixel value to be set, default: 0.0')
+    parser.add_argument('--operation', type=str, default='AND', help='How to combine multiple regions with AND or OR, default: AND')
+        
+    args = parser.parse_args()
+
+    main(args.input_image, args.region, outfile = args.outfile, inverse = args.inverse, blankval = args.blankval, operation = args.operation)
\ No newline at end of file
diff --git a/scripts/download_skymodel_target.py b/scripts/download_skymodel_target.py
index 48578e8f..4c4d15b2 100755
--- a/scripts/download_skymodel_target.py
+++ b/scripts/download_skymodel_target.py
@@ -59,7 +59,7 @@ def input2strlist_nomapfile(invar):
 
 
 ########################################################################
-def main(ms_input, SkymodelPath, Radius="5.", DoDownload="True", Source="TGSS", targetname = "pointing"):
+def main(ms_input, SkymodelPath, Radius="5.", DoDownload="True", Source="TGSS", targetname = "pointing", fluxlimit = None):
     """
     Download the skymodel for the target field
 
@@ -128,6 +128,8 @@ def main(ms_input, SkymodelPath, Radius="5.", DoDownload="True", Source="TGSS",
 
     # Treat all sources as one group (direction)
     skymodel = lsmtool.load(SkymodelPath)
+    if fluxlimit:
+        skymodel.remove('I<' + str(fluxlimit))
     skymodel.group('single', root = targetname)
     skymodel.write(clobber=True)
     
@@ -151,10 +153,12 @@ if __name__ == '__main__':
                         help='Download or not the TGSS skymodel or GSM ("Force" or "True" or "False").')
     parser.add_argument('--targetname', type=str, default='pointing',
                         help='Name of the patch of the skymodel')
+    parser.add_argument('--fluxlimit', type=float, default=None,
+                        help='Remove sources from the skymodel below the given fluxlimit [Jy],')
 
     args = parser.parse_args()
     radius=5
     if args.Radius:
         radius=args.Radius
 
-    main(args.MSfile, args.SkyTar, str(radius), args.DoDownload, args.Source, args.targetname)
+    main(args.MSfile, args.SkyTar, str(radius), args.DoDownload, args.Source, args.targetname, args.fluxlimit)
diff --git a/scripts/make_beam_reg.py b/scripts/make_beam_reg.py
new file mode 100755
index 00000000..b6d115e3
--- /dev/null
+++ b/scripts/make_beam_reg.py
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+
+# blank areas of a fits image given a region file
+
+# @author: Francesco de Gasperin, modified by Alexander Drabent
+
+def input2strlist_nomapfile(invar):
+   """ 
+   from bin/download_IONEX.py
+   give the list of MSs from the list provided as a string
+   """
+   str_list = None
+   if type(invar) is str:
+       if invar.startswith('[') and invar.endswith(']'):
+           str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
+       else:
+           str_list = [invar.strip(' \'\"')]
+   elif type(invar) is list:
+       str_list = [str(f).strip(' \'\"') for f in invar]
+   else:
+       raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
+   return str_list
+
+def getPhaseCentre(filename):
+    """
+    Get the phase centre (in degrees) of the first source of an MS.
+    """
+
+    from casacore import tables
+    import numpy as np
+
+    field_no = 0
+    ant_no   = 0
+    with tables.table(filename + "/FIELD", ack = False) as field_table:
+        direction = field_table.getcol("PHASE_DIR")
+    RA        = direction[ant_no, field_no, 0]
+    Dec       = direction[ant_no, field_no, 1]
+
+    if (RA < 0):
+        RA += 2 * np.pi
+
+    return (np.degrees(RA), np.degrees(Dec))
+
+def main(msinput, outfile, fwhm, pbcut = None, tonull = False):
+    """
+    Create a ds9 region of the beam
+    outfile : str
+        output file
+    pbcut : float, optional
+        diameter of the beam
+    tonull : bool, optional
+        arrive to the first null, not the FWHM
+    freq: min,max,med 
+        which frequency to use to estimate the beam size
+    """
+
+    import pyregion
+    from pyregion.parser_helper import Shape
+
+    print('Making PB region: ' + outfile)
+
+    msname = input2strlist_nomapfile(msinput)[0]
+    ra, dec = getPhaseCentre(msname)
+
+    if pbcut is None:
+        radius = fwhm / 2.
+    else:
+        radius = pbcut / 2.
+
+    if tonull: radius *= 2 # rough estimation
+
+    s = Shape('circle', None)
+    s.coord_format = 'fk5'
+    s.coord_list = [ ra, dec, radius ] # ra, dec, radius
+    s.coord_format = 'fk5'
+    s.attr = ([], {'width': '2', 'point': 'cross',
+                   'font': '"helvetica 16 normal roman"'})
+    s.comment = 'color=red text="beam"'
+
+    regions = pyregion.ShapeList([s])
+    regions.write(outfile)
+
+if __name__=='__main__':
+    import argparse
+    parser = argparse.ArgumentParser(description='Set to "blankval" all the pixels inside the given region')
+    
+    parser.add_argument('msinput', type=str, nargs='+', help='Input Measurement Set')
+    parser.add_argument('--outfile', type=str, help='Name of output file')
+    parser.add_argument('--fwhm', type=float, help="The FHWM of the image")
+    parser.add_argument('--pbcut', type=float, default=None, help='Diameter of the beam, default: None')
+    parser.add_argument('--tonull', type=bool, default=False, help='arrive to the first null, not the FWHM, default: False')
+        
+    args = parser.parse_args()
+
+    main(args.msinput, args.outfile, fwhm = args.fwhm, pbcut = args.pbcut, tonull = args.tonull)
\ No newline at end of file
diff --git a/scripts/make_mask.py b/scripts/make_mask.py
new file mode 100755
index 00000000..ed609a06
--- /dev/null
+++ b/scripts/make_mask.py
@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+
+# create a mask using bdsm of an image
+
+# @author: Francesco de Gasperin, modified by Alexander Drabent
+
+def make_mask(image_name, mask_name=None, threshpix=5, atrous_do=False, rmsbox=(100,10), adaptive_thresh=50,
+              write_srl=False, write_gaul=False, write_ds9=False, mask_combine=None, frequency = 54e6):
+
+    import sys, os
+    import numpy as np
+    from astropy.io import fits as pyfits
+    import bdsf
+
+    # wavelets are required to fit gaussians
+    if atrous_do or write_srl or write_ds9 or write_gaul: stop_at = None
+    else: stop_at = 'isl'
+
+    # DO THE SOURCE DETECTION
+    img = bdsf.process_image(image_name, rms_box=rmsbox, frequency=frequency,
+        thresh_isl=float(threshpix*4/5), thresh_pix=float(threshpix), rms_map=True, mean_map='zero', atrous_do=atrous_do, atrous_jmax=4,
+        adaptive_rms_box=True, adaptive_thresh=adaptive_thresh, rms_box_bright=(30,5),
+        flagging_opts=True, flag_maxsize_fwhm=0.5, stop_at=stop_at, quiet=True, debug=False)
+
+    # WRITE THE MASK FITS
+    if mask_name == None: mask_name = image_name+'.newmask'
+    if os.path.exists(mask_name): os.system('rm -r ' + mask_name)
+    img.export_image(img_type='island_mask', img_format='fits', outfile=mask_name, clobber=True)
+
+    # WRITE CATALOGUE
+    if write_srl:
+        img.write_catalog(format='fits', catalog_type='srl', outfile=mask_name.replace('fits','cat.fits'), clobber=True)
+    if write_gaul:
+        img.write_catalog(format='bbs', catalog_type='gaul', outfile=mask_name.replace('fits','skymodel'), clobber=True)
+    if write_ds9:
+        img.write_catalog(format='ds9', catalog_type='srl', outfile=mask_name.replace('fits','reg'), clobber=True)
+    del img
+
+    # do an pixel-by-pixel "OR" operation with a given mask
+    if not mask_combine is None:
+        print("Doing a pix-by-pix OR with %s." % mask_combine)
+        with pyfits.open(mask_combine) as fits:
+            data_comb = fits[0].data
+        with pyfits.open(mask_name) as fits:
+            data = fits[0].data
+            assert data.shape == data_comb.shape
+            data[(data_comb == 1.)] = 1.
+            fits[0].data = data
+            fits.writeto(mask_name, overwrite=True)
+
+    return mask_name
+
+if __name__=='__main__':
+    import optparse
+    opt = optparse.OptionParser(usage='%prog [-v|-V] imagename \n Francesco de Gasperin', version='1.0')
+    opt.add_option('-m', '--newmask', help='Mask name (default=imagename with mask in place of image)', default=None)
+    opt.add_option('-p', '--threshpix', help='Threshold pixel (default=5)', type='int', default=5)
+    opt.add_option('-a', '--adaptive_thresh', help='Adaprtive threshold (default=50)', type='int', default=50)
+    opt.add_option('-t', '--atrous_do', help='BDSM extended source detection (default=False)', action='store_true', default=False)
+    opt.add_option('-r', '--rmsbox', help='rms box size (default=100,10)', default='100,10')
+    opt.add_option('-s', '--write_srl', help='Write SRL skymodel (default=False)', action='store_true', default=False)
+    opt.add_option('-g', '--write_gaul', help='Write bbs gaul skymodel (default=False)', action='store_true', default=False)
+    opt.add_option('-d', '--write_ds9', help='Write ds9 regions (default=False)', action='store_true', default=False)
+    opt.add_option('-c', '--combinemask', help='Mask name of a mask to add to the found one (default=None)', default=None)
+    opt.add_option('-f', '--frequency', help='Provide central frequency of the observation in Hz (default=54e6)', default=54e6)
+    (options, args) = opt.parse_args()
+    
+    rmsbox = (int(options.rmsbox.split(',')[0]),int(options.rmsbox.split(',')[1]))
+    make_mask(args[0].rstrip('/'), options.newmask, options.threshpix, options.atrous_do, rmsbox, options.adaptive_thresh,
+              options.write_srl, options.write_gaul, options.write_ds9, options.combinemask, options.frequency)
diff --git a/scripts/sort_times_into_freqGroups.py b/scripts/sort_times_into_freqGroups.py
index 829792b4..cd4a358f 100755
--- a/scripts/sort_times_into_freqGroups.py
+++ b/scripts/sort_times_into_freqGroups.py
@@ -252,7 +252,7 @@ def main(MSfile, numSB=10, DP3fill=True, stepname=None, mergeLastGroup=False, tr
 
     nr_of_groups = len(groupnames)
     total_bandwidth = nr_of_groups * groupBW
-    results = {'filenames': filenames, 'groupnames': groupnames, 'total_bandwidth': total_bandwidth}
+    results = {'filenames': filenames, 'groupnames': groupnames, 'total_bandwidth': total_bandwidth, 'maxfreq': maxfreq, 'minfreq': minfreq}
 
     return(results)
     
diff --git a/steps/applycal.cwl b/steps/applycal.cwl
index 7d3ea394..4158aa11 100644
--- a/steps/applycal.cwl
+++ b/steps/applycal.cwl
@@ -82,6 +82,12 @@ inputs:
       position: 0
       prefix: applycal.missingantennabehavior=
       separate: false
+  - id: corrupt
+    type: boolean?
+    inputBinding:
+      position: 0
+      prefix: applycal.invert=False
+      separate: false
 outputs:
   - id: msout
     doc: Output Measurement Set
diff --git a/steps/blank_image_reg.cwl b/steps/blank_image_reg.cwl
new file mode 100644
index 00000000..5be18d9b
--- /dev/null
+++ b/steps/blank_image_reg.cwl
@@ -0,0 +1,70 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: blank_image_reg
+label: blank_image_reg
+baseCommand:
+  - blank_image_reg.py
+inputs:
+  - id: input_image
+    type: File
+    inputBinding:
+      position: 0
+    doc: Input image
+  - id: region
+    type: File?
+    inputBinding:
+      position: 1
+    doc: Region file
+  - id: outfile
+    type: string?
+    doc: Output file name
+    inputBinding:
+      prefix: '--outfile'
+      position: 2
+  - id: inverse
+    type: int?
+    inputBinding:
+      position: 2
+      prefix: '--inverse'
+    doc: reverse final combined mask
+  - id: blankval
+    type: float?
+    default: 0.0
+    inputBinding:
+      position: 2
+      prefix: '--blankval'
+    doc: pixel value to be set
+  - id: operation
+    type: string?
+    default: 'AND'
+    inputBinding:
+      position: 2
+      prefix: '--operation'
+    doc: how to combine multiple regions with AND or OR
+outputs:
+  - id: image
+    doc: input image as output
+    type: File
+    outputBinding:
+      glob: $(inputs.input_image.basename)
+  - id: output_image
+    doc: output image as given as output file name
+    type: File
+    outputBinding:
+      glob: $(inputs.outfile)
+  - id: logfile
+    type: File[]
+    outputBinding:
+      glob: 'blank_image_reg*.log'
+hints:
+ - class: DockerRequirement
+   dockerPull: astronrd/linc
+requirements:
+ - class: InplaceUpdateRequirement
+   inplaceUpdate: true 
+ - class: InitialWorkDirRequirement
+   listing:
+     - entry: $(inputs.input_image)
+       writable: true
+stderr: blank_image_reg.log
+stdout: blank_image_reg_err.log
diff --git a/steps/find_skymodel_target.cwl b/steps/find_skymodel_target.cwl
index c909cc42..309e6ac6 100644
--- a/steps/find_skymodel_target.cwl
+++ b/steps/find_skymodel_target.cwl
@@ -28,6 +28,8 @@ inputs:
     - id: targetname
       type: string?
       default: 'pointing'
+    - id: fluxlimit
+      type: float?
 requirements:
   - class: InlineJavascriptRequirement
   - class: NetworkAccess
@@ -40,12 +42,15 @@ requirements:
           import shutil
           import os
           import json
+          
+          null = None
 
           from download_skymodel_target import main as download_skymodel_target
 
           mss = sys.argv[1:]
           inputs = json.loads(r"""$(inputs)""")
           targetname = inputs['targetname']
+          fluxlimit = inputs['fluxlimit']
 
           SkymodelPath = inputs['SkymodelPath']
           if SkymodelPath is None:
@@ -59,7 +64,7 @@ requirements:
           Source = inputs['Source']
           DoDownload = str(inputs['DoDownload'])
 
-          output = download_skymodel_target(mss, SkymodelPath, Radius, DoDownload, Source, targetname)
+          output = download_skymodel_target(mss, SkymodelPath, Radius, DoDownload, Source, targetname, fluxlimit)
 
 outputs:
   - id: skymodel
@@ -75,4 +80,4 @@ hints:
   - class: DockerRequirement
     dockerPull: astronrd/linc
 
-stdout: find_skymodel_target.log
+stdout: find_skymodel_target.log
\ No newline at end of file
diff --git a/steps/getFWHM.cwl b/steps/getFWHM.cwl
new file mode 100644
index 00000000..c385f345
--- /dev/null
+++ b/steps/getFWHM.cwl
@@ -0,0 +1,71 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: getFWHM
+baseCommand:
+  - python3
+arguments:
+    - position: 0
+      valueFrom: getFWHM.py
+inputs:
+  - id: beamfreq
+    type: float?
+    default: 50e6
+  - id: msin
+    type: Directory[]
+    inputBinding:
+      position: 1
+outputs:
+  - id: fwhm
+    type: float
+    outputBinding:
+      loadContents: true
+      glob: 'out.json'
+      outputEval: '$(JSON.parse(self[0].contents).FWHM)'
+  - id: msout
+    doc: Output MS
+    type: Directory[]
+    outputBinding:
+      outputEval: $(inputs.msin)
+label: getFWHM
+
+requirements:
+  - class: InlineJavascriptRequirement
+  - class: InitialWorkDirRequirement
+    listing:
+     - entryname: getFWHM.py
+       entry: |
+        import sys
+        import json
+        import numpy
+        from casacore import tables
+
+        inputs = json.loads(r"""$(inputs)""")
+        beamfreq = inputs['beamfreq']
+        ms_list = sys.argv[1:]
+
+        # Following numbers are based at 60 MHz (old.astron.nl/radio-observatory/astronomers/lofar-imaging-capabilities-sensitivity/lofar-imaging-capabilities/lofa)
+
+        with tables.table(ms_list[0]+'/OBSERVATION', ack = False) as t:
+            antenna_set = t.getcell("LOFAR_ANTENNA_SET",0)
+        scale = 60e6/beamfreq
+
+        if 'HBA' in antenna_set:
+             fwhm = 10.0
+        elif 'OUTER' in antenna_set:
+             fwhm = 3.88 * scale
+        elif 'SPARSE' in antenna_set:
+             fwhm = 4.85 * scale
+        elif 'INNER' in antenna_set:
+             fwhm = 9.77 * scale
+
+        cwl_output = {"FWHM": fwhm}
+
+        with open('./out.json', 'w') as fp:
+             json.dump(cwl_output, fp)
+                
+hints:
+  - class: DockerRequirement
+    dockerPull: astronrd/linc
+
+stdout: getFWHM.log
+stderr: getFWHM_err.log
diff --git a/steps/make_beam_reg.cwl b/steps/make_beam_reg.cwl
new file mode 100644
index 00000000..a4e1bb7c
--- /dev/null
+++ b/steps/make_beam_reg.cwl
@@ -0,0 +1,55 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: make_beam_reg
+label: make_beam_reg
+baseCommand:
+  - make_beam_reg.py
+inputs:
+  - id: input_ms
+    type:
+      - Directory
+      - type: array
+        items: Directory
+    inputBinding:
+      position: 0
+    doc: Input image
+  - id: outfile
+    type: string
+    doc: Output region file name
+    inputBinding:
+      position: 1
+      prefix: --outfile
+  - id: fwhm
+    type: float
+    inputBinding:
+      position: 2
+      prefix: '--fwhm'
+    doc: The FHWM of the image
+  - id: pb_cut
+    type: float?
+    inputBinding:
+      position: 2
+      prefix: '--pbcut'
+    doc: Diameter of the beam
+  - id: to_null
+    type: boolean?
+    default: False
+    inputBinding:
+      position: 2
+      prefix: '--tonull'
+    doc: arrive to the first null, not the FWHM
+outputs:
+  - id: region
+    doc: output region file
+    type: File
+    outputBinding:
+      glob: $(inputs.outfile)
+  - id: logfile
+    type: File[]
+    outputBinding:
+      glob: 'make_beam_reg*.log'
+hints:
+ - class: DockerRequirement
+   dockerPull: astronrd/linc
+stderr: make_beam_reg.log
+stdout: make_beam_reg_err.log
diff --git a/steps/make_mask.cwl b/steps/make_mask.cwl
new file mode 100644
index 00000000..a9c9ff12
--- /dev/null
+++ b/steps/make_mask.cwl
@@ -0,0 +1,98 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: make_mask
+label: make_mask
+baseCommand:
+  - make_mask.py
+inputs:
+  - id: image
+    type: File
+    inputBinding:
+      position: 0
+    doc: Input image
+  - id: mask_name
+    type: string?
+    inputBinding:
+      position: 1
+      prefix: '-m'
+    doc: Mask name (default=imagename with mask in place of image)
+  - id: threshpix
+    type: int?
+    default: 5
+    doc: Threshold pixel
+    inputBinding:
+      prefix: '-p'
+      position: 1
+  - id: adaptive_thresh
+    type: int?
+    default: 50
+    inputBinding:
+      position: 1
+      prefix: '-a'
+    doc: Adaprtive threshold
+  - id: atrous_do
+    type: boolean?
+    default: false
+    inputBinding:
+      position: 1
+      prefix: '-t'
+    doc: BDSF extended source detection
+  - id: rmsbox
+    type: string?
+    default: '100,10'
+    inputBinding:
+      position: 0
+      prefix: '-r'
+    doc: rms box size
+  - id: write_srl
+    type: boolean?
+    default: false
+    inputBinding:
+      position: 1
+      prefix: '-s'
+    doc: Write SRL skymodel
+  - id: write_gaul
+    type: boolean?
+    default: false
+    inputBinding:
+      position: 1
+      prefix: '-g'
+    doc: Write BBS gaul skymodel
+  - id: write_ds9
+    type: boolean?
+    default: false
+    inputBinding:
+      position: 1
+      prefix: '-d'
+    doc: Write ds9 regions
+  - id: combinemask
+    type: File?
+    inputBinding:
+      position: 1
+      prefix: '-c'
+    doc: Mask to add to the found one
+  - id: frequency
+    type: float?
+    inputBinding:
+      position: 1
+      prefix: '-f'
+outputs:
+  - id: mask
+    doc: output mask
+    type: File
+    outputBinding:
+      glob: $(inputs.mask_name)
+  - id: logfile
+    type: File[]
+    outputBinding:
+      glob: 'make_mask*.log'
+hints:
+ - class: DockerRequirement
+   dockerPull: astronrd/linc
+requirements:
+  - class: InitialWorkDirRequirement
+    listing:
+      - entry: $(inputs.image)
+        writable: false
+stderr: make_mask.log
+stdout: make_mask_err.log
diff --git a/steps/sort_times_into_freqGroups.cwl b/steps/sort_times_into_freqGroups.cwl
index 0ea9d57c..5f884603 100644
--- a/steps/sort_times_into_freqGroups.cwl
+++ b/steps/sort_times_into_freqGroups.cwl
@@ -59,10 +59,14 @@ requirements:
           filenames  = output['filenames']
           groupnames = output['groupnames']
           total_bandwidth = output['total_bandwidth']
+          midfreq    = (output['maxfreq'] + output['minfreq']) / 2.
+          minfreq    = output['minfreq']
 
           cwl_output = {}
           cwl_output['groupnames'] = groupnames
           cwl_output['total_bandwidth'] = total_bandwidth
+          cwl_output['midfreq'] = midfreq
+          cwl_output['minfreq'] = minfreq
 
           with open('./filenames.json', 'w') as fp:
               json.dump(filenames, fp)
@@ -87,6 +91,18 @@ outputs:
         loadContents: true
         glob: 'out.json'
         outputEval: $(JSON.parse(self[0].contents).total_bandwidth)
+  - id: midfreq
+    type: float
+    outputBinding:
+        loadContents: true
+        glob: 'out.json'
+        outputEval: $(JSON.parse(self[0].contents).midfreq)
+  - id: minfreq
+    type: float
+    outputBinding:
+        loadContents: true
+        glob: 'out.json'
+        outputEval: $(JSON.parse(self[0].contents).minfreq)
   - id: logfile
     type: File
     outputBinding:
diff --git a/steps/wsclean.cwl b/steps/wsclean.cwl
index 8a845542..e435d8bb 100644
--- a/steps/wsclean.cwl
+++ b/steps/wsclean.cwl
@@ -15,7 +15,8 @@ inputs:
     default:
       - 2500
       - 2500
-    type: int[]
+    type: 
+      - int[]
     inputBinding:
       position: 1
       shellQuote: false
@@ -47,11 +48,20 @@ inputs:
     default: false
     type: 
       - boolean?
-      - int?
+      - float?
     inputBinding:
       shellQuote: false
       position: 1
       prefix: '-auto-threshold'
+  - id: auto_mask
+    default: false
+    type: 
+      - boolean?
+      - float?
+    inputBinding:
+      shellQuote: false
+      position: 1
+      prefix: '-auto-mask'
   - id: multiscale
     default: false
     type: 
@@ -60,6 +70,22 @@ inputs:
       position: 1
       shellQuote: false
       prefix: '-multiscale'
+  - id: multiscale_scales
+    default: false
+    type: 
+      - boolean?
+      - string?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-multiscale-scales'
+  - id: multiscale_scale_bias
+    default: false
+    type: float?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-multiscale-scale-bias'
   - id: mgain
     default: false
     type: 
@@ -173,6 +199,24 @@ inputs:
       position: 1
       prefix: '-maxuvw-m'
       shellQuote: false
+  - id: minuv-l
+    default: false
+    type:
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      prefix: '-minuv-l'
+      shellQuote: false
+  - id: maxuv-l
+    default: false
+    type:
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      prefix: '-maxuv-l'
+      shellQuote: false
   - id: tempdir
     type: 
       - boolean?
@@ -182,13 +226,43 @@ inputs:
       position: 1
       prefix: '-temp-dir'
       shellQuote: false
-  - id: model_update
-    default: true
+  - id: no_model_update
+    default: false
     type: boolean?
     inputBinding:
       position: 1
       prefix: '-no-update-model-required'
       shellQuote: false
+  - id: no_fit_beam
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-no-fit-beam'
+      shellQuote: false
+  - id: circular_beam
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-circular-beam'
+      shellQuote: false
+  - id: beam_size
+    default: false
+    type:
+      - boolean?
+      - float?
+    inputBinding:
+      position: 1
+      prefix: '-beam-size'
+      shellQuote: false
+  - id: local_rms
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-local-rms'
+      shellQuote: false
   - id: image_name
     default: 'pointing'
     type: string?
@@ -196,31 +270,97 @@ inputs:
       position: 1
       prefix: '-name'
       shellQuote: false
+  - id: save_source_list
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-save-source-list'
+      shellQuote: false
+  - id: do_predict
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-predict'
+      shellQuote: false
+  - id: baseline_averaging
+    default: false
+    type: float?
+    inputBinding:
+      position: 1
+      prefix: -baseline-averaging
+  - id: reuse_dirty
+    type: File?
+    inputBinding:
+      position: 1
+      prefix: '-reuse-dirty'
+      shellQuote: false
+      valueFrom: $(self.nameroot.replaceAll('-dirty',''))
+  - id: reuse_psf
+    type: File?
+    inputBinding:
+      position: 1
+      prefix: '-reuse-psf'
+      shellQuote: false
+      valueFrom: $(self.nameroot.replaceAll('-psf',''))
+  - id: fits_mask
+    type: File?
+    inputBinding:
+      position: 1
+      prefix: '-fits-mask'
+      shellQuote: false
+  - id: fits_model
+    type: File?
+  - id: fits_image
+    type: File?
 outputs:
   - id: dirty_image
-    type: File
+    type: File?
     outputBinding:
       glob: [$(inputs.image_name)-MFS-dirty.fits, $(inputs.image_name)-dirty.fits]
+  - id: psf_image
+    type: File?
+    outputBinding:
+      glob: [$(inputs.image_name)-MFS-psf.fits, $(inputs.image_name)-psf.fits]
   - id: image
     type: File
     outputBinding:
       glob: [$(inputs.image_name)-MFS-image.fits, $(inputs.image_name)-image.fits]
+  - id: model
+    type: File?
+    outputBinding:
+      glob: [$(inputs.image_name)-MFS-model.fits, $(inputs.image_name)-model.fits]
   - id: logfile
     type: File[]
     outputBinding:
       glob: 'wsclean*.log'
+  - id: msout
+    type: Directory[]
+    outputBinding:
+      outputEval: $(inputs.msin)
 label: WSClean
 hints:
   - class: DockerRequirement
     dockerPull: 'astronrd/linc'
+requirements:
+  - class: InplaceUpdateRequirement
+    inplaceUpdate: true 
   - class: InitialWorkDirRequirement
     listing:
-      - $(inputs.msin)
-requirements:
+      - entry: $(inputs.msin)
+        writable: true
+      - entry: $(inputs.reuse_psf)
+        writable: false
+      - entry: $(inputs.reuse_dirty)
+        writable: false
+      - entry: $(inputs.fits_model)
+        writable: false
+      - entry: $(inputs.fits_image)
+        writable: false
   - class: ShellCommandRequirement
   - class: InlineJavascriptRequirement
     expressionLib:
       - { $include: 'utils.js' }
-
 stdout: wsclean.log
 stderr: wsclean_err.log
\ No newline at end of file
diff --git a/steps/wsclean_predict.cwl b/steps/wsclean_predict.cwl
new file mode 100644
index 00000000..1787abfb
--- /dev/null
+++ b/steps/wsclean_predict.cwl
@@ -0,0 +1,360 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: wsclean_predict
+baseCommand:
+  - wsclean
+inputs:
+  - id: msin
+    type: Directory[]
+    inputBinding:
+      position: 2
+      shellQuote: false
+      itemSeparator: ','
+      valueFrom: $(concatenate_path_wsclean(self))
+  - id: image_size
+    default:
+      - 2500
+      - 2500
+    type: 
+      - int[]
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-size'
+  - id: image_scale
+    default: '36asec'
+    type: string?
+    inputBinding:
+      position: 1
+      prefix: '-scale'
+      shellQuote: false
+  - id: niter
+    default: 1000
+    type: int?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-niter'
+  - id: nmiter
+    default: false
+    type: 
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-nmiter'
+  - id: auto_threshold
+    default: false
+    type: 
+      - boolean?
+      - float?
+    inputBinding:
+      shellQuote: false
+      position: 1
+      prefix: '-auto-threshold'
+  - id: auto_mask
+    default: false
+    type: 
+      - boolean?
+      - float?
+    inputBinding:
+      shellQuote: false
+      position: 1
+      prefix: '-auto-mask'
+  - id: multiscale
+    default: false
+    type: 
+      - boolean?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-multiscale'
+  - id: multiscale_scales
+    default: false
+    type: 
+      - boolean?
+      - string?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-multiscale-scales'
+  - id: multiscale_scale_bias
+    default: false
+    type: float?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-multiscale-scale-bias'
+  - id: mgain
+    default: false
+    type: 
+      - boolean?
+      - float?
+    inputBinding:
+      shellQuote: false
+      position: 1
+      prefix: '-mgain'
+  - id: ncpu
+    default: false
+    type: 
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-j'
+  - id: parallel-gridding
+    default: false
+    type: 
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-parallel-gridding'
+  - id: parallel-deconvolution
+    default: false
+    type: 
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-parallel-deconvolution'
+  - id: parallel-reordering
+    default: false
+    type: 
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-parallel-reordering'
+  - id: channels-out
+    default: false
+    type: 
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-channels-out'
+  - id: deconvolution-channels
+    default: false
+    type: 
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-deconvolution-channels'
+  - id:  fit-spectral-pol
+    default: false
+    type: 
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-fit-spectral-pol'
+  - id: join-channels
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-join-channels'
+  - id: use-wgridder
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-use-wgridder'
+  - id: taper-gaussian
+    default: false
+    type: 
+      - boolean?
+      - string?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-taper-gaussian'
+  - id: weighting
+    default: false
+    type: 
+      - boolean?
+      - string?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-weight'
+  - id: maxuvw-m
+    default: false
+    type: 
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      prefix: '-maxuvw-m'
+      shellQuote: false
+  - id: minuv-l
+    default: false
+    type:
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      prefix: '-minuv-l'
+      shellQuote: false
+  - id: maxuv-l
+    default: false
+    type:
+      - boolean?
+      - int?
+    inputBinding:
+      position: 1
+      prefix: '-maxuv-l'
+      shellQuote: false
+  - id: tempdir
+    type: 
+      - boolean?
+      - string?
+    default: false
+    inputBinding:
+      position: 1
+      prefix: '-temp-dir'
+      shellQuote: false
+  - id: no_model_update
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-no-update-model-required'
+      shellQuote: false
+  - id: no_fit_beam
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-no-fit-beam'
+      shellQuote: false
+  - id: circular_beam
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-circular-beam'
+      shellQuote: false
+  - id: beam_size
+    default: false
+    type:
+      - boolean?
+      - float?
+    inputBinding:
+      position: 1
+      prefix: '-beam-size'
+      shellQuote: false
+  - id: local_rms
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-local-rms'
+      shellQuote: false
+  - id: image_name
+    default: 'pointing'
+    type: string?
+    inputBinding:
+      position: 1
+      prefix: '-name'
+      shellQuote: false
+  - id: save_source_list
+    default: false
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-save-source-list'
+      shellQuote: false
+  - id: do_predict
+    default: true
+    type: boolean?
+    inputBinding:
+      position: 1
+      prefix: '-predict'
+      shellQuote: false
+  - id: baseline_averaging
+    default: false
+    type: float?
+    inputBinding:
+      position: 1
+      prefix: -baseline-averaging
+  - id: reuse_dirty
+    type: File?
+    inputBinding:
+      position: 1
+      prefix: '-reuse-dirty'
+      shellQuote: false
+      valueFrom: $(self.nameroot.replaceAll('-dirty',''))
+  - id: reuse_psf
+    type: File?
+    inputBinding:
+      position: 1
+      prefix: '-reuse-psf'
+      shellQuote: false
+      valueFrom: $(self.nameroot.replaceAll('-psf',''))
+  - id: fits_mask
+    type: File?
+    inputBinding:
+      position: 1
+      prefix: '-fits-mask'
+      shellQuote: false
+  - id: fits_model
+    type: File?
+  - id: fits_image
+    type: File?
+outputs:
+  - id: dirty_image
+    type: File?
+    outputBinding:
+      glob: [$(inputs.image_name)-MFS-dirty.fits, $(inputs.image_name)-dirty.fits]
+  - id: psf_image
+    type: File?
+    outputBinding:
+      glob: [$(inputs.image_name)-MFS-psf.fits, $(inputs.image_name)-psf.fits]
+  - id: model
+    type: File?
+    outputBinding:
+      glob: [$(inputs.image_name)-MFS-model.fits, $(inputs.image_name)-model.fits]
+  - id: logfile
+    type: File[]
+    outputBinding:
+      glob: 'wsclean*.log'
+  - id: msout
+    type: Directory[]
+    outputBinding:
+      outputEval: $(inputs.msin)
+label: WSClean
+hints:
+  - class: DockerRequirement
+    dockerPull: 'astronrd/linc'
+requirements:
+  - class: InitialWorkDirRequirement
+    listing:
+      - entry: $(inputs.msin)
+        writable: true
+      - entry: $(inputs.reuse_psf)
+        writable: false
+      - entry: $(inputs.reuse_dirty)
+        writable: false
+      - entry: $(inputs.fits_model)
+        writable: false
+      - entry: $(inputs.fits_image)
+        writable: false
+  - class: ShellCommandRequirement
+  - class: InlineJavascriptRequirement
+    expressionLib:
+      - { $include: 'utils.js' }
+stdout: wsclean_predict.log
+stderr: wsclean_predict_err.log
\ No newline at end of file
diff --git a/workflows/HBA_target.cwl b/workflows/HBA_target.cwl
index dcd921fb..d7969eb8 100644
--- a/workflows/HBA_target.cwl
+++ b/workflows/HBA_target.cwl
@@ -165,6 +165,10 @@ inputs:
     default: 0.0
   - id: wsclean_tmpdir
     type: string?
+  - id: selfcal_region
+    type: File?
+  - id: skymodel_fluxlimit
+    type: float?
 outputs:
   - id: calibrated_data
     outputSource:
@@ -297,6 +301,10 @@ steps:
         source: chunkduration
       - id: wsclean_tmpdir
         source: wsclean_tmpdir
+      - id: selfcal_region
+        source: selfcal_region
+      - id: skymodel_fluxlimit
+        source: skymodel_fluxlimit
     out:
       - id: logfiles
       - id: msout
diff --git a/workflows/LBA_target.cwl b/workflows/LBA_target.cwl
index b01f0d25..953f287e 100644
--- a/workflows/LBA_target.cwl
+++ b/workflows/LBA_target.cwl
@@ -166,6 +166,11 @@ inputs:
     default: 3600.0
   - id: wsclean_tmpdir
     type: string?
+  - id: selfcal_region
+    type: File?
+  - id: skymodel_fluxlimit
+    type: float?
+    default: 1.0
 outputs:
   - id: calibrated_data
     outputSource:
@@ -298,6 +303,10 @@ steps:
         source: chunkduration
       - id: wsclean_tmpdir
         source: wsclean_tmpdir
+      - id: selfcal_region
+        source: selfcal_region
+      - id: skymodel_fluxlimit
+        source: skymodel_fluxlimit
     out:
       - id: logfiles
       - id: msout
diff --git a/workflows/linc_target.cwl b/workflows/linc_target.cwl
index 594da77c..9223419f 100644
--- a/workflows/linc_target.cwl
+++ b/workflows/linc_target.cwl
@@ -165,6 +165,10 @@ inputs:
     default: 0.0
   - id: wsclean_tmpdir
     type: string?
+  - id: selfcal_region
+    type: File?
+  - id: skymodel_fluxlimit
+    type: float?
 outputs:
   - id: inspection
     outputSource:
@@ -216,6 +220,10 @@ steps:
           - flag_baselines
       - id: process_baselines_target
         source: process_baselines_target
+      - id: num_SBs_per_group
+        source: num_SBs_per_group
+      - id: reference_stationSB
+        source: reference_stationSB
       - id: filter_baselines
         source: filter_baselines
       - id: raw_data
@@ -283,6 +291,8 @@ steps:
         source: lbfgs_historysize
       - id: lbfgs_robustdof
         source: lbfgs_robustdof
+      - id: skymodel_fluxlimit
+        source: skymodel_fluxlimit
     out:
       - id: compare_stations_filter
       - id: outh5parm
@@ -293,6 +303,10 @@ steps:
       - id: initial_flags_join_out
       - id: logfiles
       - id: target_sourcedb
+      - id: total_bandwidth
+      - id: midfreq
+      - id: groupnames
+      - id: filenames
     run: ./linc_target/prep.cwl
     label: prep
   - id: gsmcal
@@ -303,12 +317,10 @@ steps:
         source: prep/msout
       - id: filter_baselines
         source: process_baselines_target
-      - id: num_SBs_per_group
-        source: num_SBs_per_group
-      - id: reference_stationSB
-        source: reference_stationSB
       - id: target_skymodel
         source: prep/target_sourcedb
+      - id: targetname
+        source: get_targetname/targetname
       - id: do_smooth
         source: do_smooth
       - id: propagatesolutions
@@ -341,6 +353,18 @@ steps:
         source: selfcal
       - id: chunkduration
         source: chunkduration
+      - id: wsclean_tmpdir
+        source: wsclean_tmpdir
+      - id: total_bandwidth
+        source: prep/total_bandwidth
+      - id: midfreq
+        source: prep/midfreq
+      - id: groupnames
+        source: prep/groupnames
+      - id: filenames
+        source: prep/filenames
+      - id: selfcal_region
+        source: selfcal_region
     out:
       - id: msout
       - id: outh5parm
@@ -351,7 +375,6 @@ steps:
       - id: inspection
       - id: logfiles
       - id: removed_bands
-      - id: total_bandwidth
       - id: out_refant
     run: ./linc_target/gsmcal.cwl
     label: gsmcal
@@ -377,7 +400,7 @@ steps:
       - id: skymodel_source
         source: skymodel_source
       - id: total_bandwidth
-        source: gsmcal/total_bandwidth
+        source: prep/total_bandwidth
       - id: check_Ateam_separation.json
         source: prep/check_Ateam_separation.json
       - id: filter_baselines
diff --git a/workflows/linc_target/corrupt_model.cwl b/workflows/linc_target/corrupt_model.cwl
new file mode 100644
index 00000000..f482a025
--- /dev/null
+++ b/workflows/linc_target/corrupt_model.cwl
@@ -0,0 +1,167 @@
+class: Workflow
+cwlVersion: v1.2
+id: corrupt_model
+label: corrupt_model
+inputs:
+  - id: msin
+    type: Directory
+  - id: max_dp3_threads
+    type: int?
+  - id: parmdb
+    type: File
+  - id: skymodel_source
+    type: string?
+    default: 'GSM'
+outputs:
+  - id: msout
+    outputSource:
+      - subtract_model/msout
+    type: Directory
+  - id: logfile
+    outputSource: concat_logfiles_corrupt/output
+    type: File
+
+steps:
+  - id: corrupt_slowtec
+    in:
+      - id: max_dp3_threads
+        source: max_dp3_threads
+      - id: msin
+        source: msin
+      - id: msin_datacolumn
+        default: MODEL_DATA
+      - id: parmdb
+        source: parmdb
+      - id: msout_datacolumn
+        default: MODEL_DATA
+      - id: storagemanager
+        default: Dysco
+      - id: databitrate
+        default: 0
+      - id: correction
+        source:
+          - skymodel_source
+        valueFrom: $("slow"+self+"tec1")
+      - id: solset
+        default: target
+      - id: missingantennabehavior
+        default: 'unit'
+      - id: corrupt
+        default: true
+    out:
+      - id: msout
+      - id: logfile
+    run: ../../steps/applycal.cwl
+    label: corrupt_slowtec
+  - id: corrupt_tec
+    in:
+      - id: max_dp3_threads
+        source: max_dp3_threads
+      - id: msin
+        source: corrupt_slowtec/msout
+      - id: msin_datacolumn
+        default: MODEL_DATA
+      - id: parmdb
+        source: parmdb
+      - id: msout_datacolumn
+        default: MODEL_DATA
+      - id: storagemanager
+        default: Dysco
+      - id: databitrate
+        default: 0
+      - id: correction
+        source:
+          - skymodel_source
+        valueFrom: $(self+"tec1")
+      - id: solset
+        default: target
+      - id: corrupt
+        default: true
+    out:
+      - id: msout
+      - id: logfile
+    run: ../../steps/applycal.cwl
+    label: corrupt_tec
+  - id: corrupt_fr
+    in:
+      - id: max_dp3_threads
+        source: max_dp3_threads
+      - id: msin
+        source: corrupt_tec/msout
+      - id: msin_datacolumn
+        default: MODEL_DATA
+      - id: parmdb
+        source: parmdb
+      - id: msout_datacolumn
+        default: MODEL_DATA
+      - id: storagemanager
+        default: Dysco
+      - id: databitrate
+        default: 0
+      - id: correction
+        default: "faraday"
+      - id: solset
+        default: target
+      - id: corrupt
+        default: true
+    out:
+      - id: msout
+      - id: logfile
+    run: ../../steps/applycal.cwl
+    label: corrupt_fr
+  - id: corrupt_amp
+    in:
+      - id: max_dp3_threads
+        source: max_dp3_threads
+      - id: msin
+        source: corrupt_fr/msout
+      - id: msin_datacolumn
+        default: MODEL_DATA
+      - id: parmdb
+        source: parmdb
+      - id: msout_datacolumn
+        default: MODEL_DATA
+      - id: storagemanager
+        default: Dysco
+      - id: databitrate
+        default: 0
+      - id: correction
+        source: skymodel_source
+        valueFrom: $(self+"amplitude")
+      - id: solset
+        default: target
+      - id: corrupt
+        default: true
+    out:
+      - id: msout
+      - id: logfile
+    run: ../../steps/applycal.cwl
+    label: corrupt_amp
+  - id: subtract_model
+    in:
+      - id: msin
+        source: corrupt_amp/msout
+      - id: command
+        default: "set CORRECTED_DATA = DATA - MODEL_DATA"
+    out:
+      - id: msout
+      - id: logfile
+    run: ../../steps/taql.cwl
+    label: subtract_model
+  - id: concat_logfiles_corrupt
+    in:
+      - id: file_list
+        source:
+          - corrupt_slowtec/logfile
+          - corrupt_tec/logfile
+          - corrupt_fr/logfile
+          - corrupt_amp/logfile
+          - subtract_model/logfile
+        linkMerge: merge_flattened
+      - id: file_prefix
+        default: corrupt_model
+    out:
+      - id: output
+    run: ../../steps/concatenate_files.cwl
+    label: concat_logfiles_corrupt
+requirements: []
diff --git a/workflows/linc_target/finalize.cwl b/workflows/linc_target/finalize.cwl
index cf6cb38c..835ccfd1 100644
--- a/workflows/linc_target/finalize.cwl
+++ b/workflows/linc_target/finalize.cwl
@@ -309,6 +309,8 @@ steps:
         source: wsclean_tmpdir
       - id: image_name
         source: targetname
+      - id: no_model_update
+        default: true
     out:
       - id: dirty_image
       - id: image
diff --git a/workflows/linc_target/gsmcal.cwl b/workflows/linc_target/gsmcal.cwl
index b9fac45f..38e944da 100644
--- a/workflows/linc_target/gsmcal.cwl
+++ b/workflows/linc_target/gsmcal.cwl
@@ -10,12 +10,6 @@ inputs:
   - id: filter_baselines
     type: string?
     default: '[CR]S*&'
-  - id: num_SBs_per_group
-    type: int?
-    default: 10
-  - id: reference_stationSB
-    type: int?
-    default: null
   - id: target_skymodel
     type: 
       - File
@@ -68,6 +62,24 @@ inputs:
   - id: chunkduration
     type: float?
     default: 0.0
+  - id: targetname
+    type: string?
+    default: 'pointing'
+  - id: wsclean_tmpdir
+    type: string?
+    default: '/tmp'
+  - id: total_bandwidth
+    type: float?
+    default: 50
+  - id: midfreq
+    type: float?
+    default: 50e6
+  - id: groupnames
+    type: string[]
+  - id: filenames
+    type: File
+  - id: selfcal_region
+    type: File?
 outputs:
   - id: msout
     outputSource:
@@ -118,7 +130,6 @@ outputs:
     outputSource:
       - concat_logfiles_identify/output
       - concat_logfiles_RefAnt/output
-      - sort_times_into_freqGroups/logfile
       - concat_logfiles_calib/output
       - concat_logfiles_dp3concat/output
       - concat_logfiles_blsmooth/output
@@ -132,10 +143,6 @@ outputs:
     outputSource:
       - check_unflagged_fraction/filenames
     type: string[]
-  - id: total_bandwidth
-    outputSource:
-      - sort_times_into_freqGroups/total_bandwidth
-    type: float
 steps:
   - id: identifybadantennas
     in:
@@ -185,28 +192,7 @@ steps:
       - id: flagged_fraction_antenna
     run: ../../steps/findRefAnt_join.cwl
     label: Ateam_flags_join
-  - id: sort_times_into_freqGroups
-    in:
-      - id: msin
-        source:
-          - msin
-      - id: numbands
-        source: num_SBs_per_group
-      - id: DP3fill
-        default: true
-      - id: stepname
-        default: .dp3concat
-      - id: firstSB
-        source: reference_stationSB
-      - id: truncateLastSBs
-        default: false
-    out:
-      - id: filenames
-      - id: groupnames
-      - id: total_bandwidth
-      - id: logfile
-    run: ../../steps/sort_times_into_freqGroups.cwl
-    label: sorttimesintofreqGroups
+
   - id: concat_logfiles_dp3concat
     in:
       - id: file_list
@@ -307,12 +293,12 @@ steps:
   - id: concat
     in:
       - id: msin
-        source:
+        source: 
           - msin
       - id: group_id
-        source: sort_times_into_freqGroups/groupnames
+        source: groupnames
       - id: groups_specification
-        source: sort_times_into_freqGroups/filenames
+        source: filenames
       - id: filter_baselines
         source: identifybadantennas_join/filter_out
       - id: avg_timeresolution_concat
@@ -446,6 +432,24 @@ steps:
         source: identifybadantennas_join/filter_out
       - id: insolutions
         source: insolutions
+      - id: targetname
+        source: targetname
+      - id: wsclean_tmpdir
+        source: wsclean_tmpdir
+      - id: total_bandwidth
+        source: total_bandwidth
+      - id: midfreq
+        source: midfreq
+      - id: selfcal_region
+        source: selfcal_region
+      - id: rfistrategy
+        source: rfistrategy
+      - id: aoflag_reorder
+        source: aoflag_reorder
+      - id: aoflag_chunksize
+        source: aoflag_chunksize
+      - id: aoflag_freqconcat
+        source: aoflag_freqconcat
       - id: execute
         source: selfcal
     out:
diff --git a/workflows/linc_target/imaging_subtract.cwl b/workflows/linc_target/imaging_subtract.cwl
new file mode 100644
index 00000000..c3e3b193
--- /dev/null
+++ b/workflows/linc_target/imaging_subtract.cwl
@@ -0,0 +1,666 @@
+class: Workflow
+cwlVersion: v1.2
+id: imaging_subtract
+label: imaging_subtract
+inputs:
+  - id: max_dp3_threads
+    type: int?
+  - id: msin
+    type: Directory[]
+  - id: targetname
+    type: string?
+    default: 'pointing'
+  - id: wsclean_tmpdir
+    type: string?
+    default: '/tmp'
+  - id: total_bandwidth
+    type: float?
+    default: 50
+  - id: midfreq
+    type: float?
+    default: 50e6
+  - id: selfcal_region
+    type: File?
+  - id: rfistrategy
+    type:
+      - File?
+      - string?
+  - id: aoflag_reorder
+    type: boolean?
+    default: false
+  - id: aoflag_chunksize
+    type: int?
+    default: 2000
+  - id: aoflag_freqconcat
+    type: boolean?
+    default: true
+  - id: parmdb
+    type: File
+  - id: skymodel_source
+    type: string?
+    default: 'GSM'
+outputs:
+  - id: msout
+    outputSource:
+      - recreate_model/msout
+    type: Directory[]
+  - id: logfile
+    outputSource:
+      - concat_logfiles_predict/output
+      - concat_logfiles_subtract/output
+      - concat_logfiles_lowres/output
+      - concat_logfiles_aoflag/output
+      - concat_logfiles_large/output
+      - concat_logfiles_corrupt/output
+      - concat_logfiles_recreate/output
+    linkMerge: merge_flattened
+    type: File[]
+    pickValue: all_non_null
+  - id: inspection
+    outputSource:
+      - make_mask_image/image
+      - make_mask/mask
+      - blank_image_reg/image
+      - blank_image_reg_1/image
+      - make_beam_reg/region
+      - blank_image_reg_2/output_image
+      - image_lowres/image
+      - image_large/image
+    linkMerge: merge_flattened
+    type: File[]
+    pickValue: all_non_null
+steps:
+  - id: getFWHM
+    in:
+      - id: msin
+        source: msin
+      - id: beamfreq
+        source: midfreq
+    out:
+      - id: fwhm
+      - id: msout
+    run: ../../steps/getFWHM.cwl
+    label: getFWHM
+  - id: make_mask_image
+    in:
+      - id: msin
+        source: getFWHM/msout # workaround for CWL bug
+      - id: image_name
+        source: targetname
+        valueFrom: '$(self+"_mask")'
+      - id: image_size
+        source: getFWHM/fwhm
+        valueFrom: '$((Math.round(2.1 * self * 3600 / 10.)) % 2 > 0 ? [Math.round(2.1 * self * 3600 / 10.) + 1 , Math.round(2.1 * self * 3600 / 10.) + 1] : [Math.round(2.1 * self * 3600 / 10.), Math.round(2.1 * self * 3600 / 10.)])'
+      - id: image_scale
+        default: '10arcsec'
+      - id: weighting
+        default: 'briggs -0.3'
+      - id: niter
+        default: 1000000
+      - id: minuv-l
+        default: 30
+      - id: parallel-gridding
+        default: 2
+      - id: maxuv-l
+        default: 4500
+      - id: mgain
+        default: 0.85
+      - id: parallel-deconvolution
+        default: 512
+      - id: local_rms
+        default: true
+      - id: auto_threshold
+        default: 4
+      - id: join-channels
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/4.e6) > 1)'
+      - id: fit-spectral-pol
+        source: total_bandwidth
+        valueFrom: '$(self > 25e6 ? 5 : Math.round(self/4.e6) > 1 ? 3 : false)'
+      - id: channels-out
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/4.e6) > 0 ? Math.round(self/4.e6) : 1)'
+      - id: deconvolution-channels
+        source: total_bandwidth
+        valueFrom: '$(self > 25e6 ? 5 : Math.round(self/4.e6) > 1 ? 3 : false)'
+      - id: use-wgridder
+        default: true
+      - id: tempdir
+        source: wsclean_tmpdir
+    out:
+      - id: dirty_image
+      - id: psf_image
+      - id: model
+      - id: image
+      - id: logfile
+      - id: msout
+    run: ../../steps/wsclean.cwl
+    label: make_mask_image
+  - id: make_mask
+    in:
+      - id: image
+        source: make_mask_image/image
+      - id: mask_name
+        source: targetname
+        valueFrom: '$(self+".mask")'
+      - id: atrous_do
+        default: true
+      - id: threshpix
+        default: 5
+      - id: frequency
+        source: midfreq
+    out:
+      - id: mask
+      - id: logfile
+    run: ../../steps/make_mask.cwl
+    label: make_mask
+  - id: blank_image_reg
+    in:
+      - id: input_image
+        source: make_mask/mask
+      - id: region
+        source: selfcal_region
+      - id: blankval
+        default: 1.0
+    out:
+      - id: image
+      - id: logfile
+    run: ../../steps/blank_image_reg.cwl
+    label: blank_image_reg
+    when: $(inputs.region != null)
+  - id: image_predict
+    in:
+      - id: msin
+        source: make_mask_image/msout # workaround for CWL bug
+      - id: image_name
+        source: targetname
+        valueFrom: '$(self+"_mask")'
+      - id: save_source_list
+        default: true
+      - id: image_size
+        source: getFWHM/fwhm
+        valueFrom: '$((Math.round(2.1 * self * 3600 / 10.)) % 2 > 0 ? [Math.round(2.1 * self * 3600 / 10.) + 1 , Math.round(2.1 * self * 3600 / 10.) + 1] : [Math.round(2.1 * self * 3600 / 10.), Math.round(2.1 * self * 3600 / 10.)])'
+      - id: image_scale
+        default: '10arcsec'
+      - id: weighting
+        default: 'briggs -0.3'
+      - id: no_model_update
+        default: true
+      - id: niter
+        default: 1000000
+      - id: minuv-l
+        default: 30
+      - id: parallel-gridding
+        default: 2
+      - id: baseline_averaging
+        default: getFWHM/fwhm
+        valueFrom: '$(1.87e3 * 60000. * 2. * Math.PI / (24 * 60 * 60 * Math.round(2.1 * self * 3600 / 10.)) < 1 ? false : 1.87e3 * 60000. * 2. * Math.PI / (24 * 60 * 60 * Math.round(2.1 * self * 3600 / 10.)) > 10 ? 10 : 1.87e3 * 60000. * 2. * Math.PI / (24 * 60 * 60 * Math.round(2.1 * self * 3600 / 10.)))'
+      - id: maxuv-l
+        default: 4500
+      - id: mgain
+        default: 0.85
+      - id: parallel-deconvolution
+        default: 512
+      - id: auto_threshold
+        default: 3
+      - id: fits_mask
+        source:
+          - make_mask/mask
+          - blank_image_reg/image
+        pickValue: first_non_null
+      - id: fits_model
+        source: make_mask_image/model
+      - id: fits_image
+        source: make_mask_image/image
+      - id: join-channels
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/4.e6) > 1)'
+      - id: fit-spectral-pol
+        source: total_bandwidth
+        valueFrom: '$(self > 25e6 ? 5 : Math.round(self/4.e6) > 1 ? 3 : false)'
+      - id: channels-out
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/4.e6) > 0 ? Math.round(self/4.e6) : 1)'
+      - id: multiscale
+        default: true
+      - id: multiscale_scale_bias
+        default: 0.6
+      - id: deconvolution-channels
+        source: total_bandwidth
+        valueFrom: '$(self > 25e6 ? 5 : Math.round(self/4.e6) > 1 ? 3 : false)'
+      - id: use-wgridder
+        default: true
+      - id: tempdir
+        source: wsclean_tmpdir
+      - id: do_predict
+        default: true
+      - id: reuse_dirty
+        source: make_mask_image/dirty_image
+      - id: reuse_psf
+        source: make_mask_image/psf_image
+    out:
+      - id: msout
+      - id: image
+      - id: model
+      - id: logfile
+    run: ../../steps/wsclean.cwl
+    label: image_predict
+  - id: subtract_model
+    in:
+      - id: msin
+        source: image_predict/msout
+      - id: command
+        default: "set CORRECTED_DATA = CORRECTED_DATA - MODEL_DATA"
+    out:
+      - id: msout
+      - id: logfile
+    run: ../../steps/taql.cwl
+    label: subtract_model
+    scatter:
+      - msin
+  - id: image_tmp
+    in:
+      - id: msin
+        source: subtract_model/msout
+      - id: image_name
+        source: targetname
+        valueFrom: '$(self+"_tmp")'
+      - id: image_size
+        source: getFWHM/fwhm
+        valueFrom: '$((Math.round(2.1 * self * 3600 / 10.)) % 2 > 0 ? [Math.round(2.1 * self * 3600 / 10.) + 1 , Math.round(2.1 * self * 3600 / 10.) + 1] : [Math.round(2.1 * self * 3600 / 10.), Math.round(2.1 * self * 3600 / 10.)])'
+      - id: image_scale
+        default: '30arcsec'
+      - id: use-wgridder
+        default: true
+      - id: tempdir
+        source: wsclean_tmpdir
+    out:
+      - id: msout
+      - id: image
+      - id: logfile
+    run: ../../steps/wsclean.cwl
+    label: image_tmp
+  - id: make_beam_reg
+    in:
+      - id: input_ms
+        source: msin
+      - id: outfile
+        source: targetname
+        valueFrom: '$(self+".reg")'
+      - id: fwhm
+        source: getFWHM/fwhm
+    out:
+      - id: region
+      - id: logfile
+    run: ../../steps/make_beam_reg.cwl
+    label: make_beam_reg
+  - id: blank_image_reg_1
+    in:
+      - id: input_image
+        source: image_tmp/image
+      - id: outfile
+        source: targetname
+        valueFrom: '$(self+"_blank1.fits")'
+      - id: region
+        source: make_beam_reg/region
+      - id: blankval
+        default: 0.0
+    out:
+      - id: image
+      - id: output_image
+      - id: logfile
+    run: ../../steps/blank_image_reg.cwl
+    label: blank_image_reg_1
+  - id: blank_image_reg_2
+    in:
+      - id: input_image
+        source: blank_image_reg_1/output_image
+      - id: outfile
+        source: targetname
+        valueFrom: '$(self+"_blank2.fits")'
+      - id: region
+        source: make_beam_reg/region
+      - id: blankval
+        default: 1.0
+      - id: inverse
+        default: 1
+    out:
+      - id: output_image
+      - id: logfile
+    run: ../../steps/blank_image_reg.cwl
+    label: blank_image_reg_2
+  - id: image_lowres
+    in:
+      - id: msin
+        source: image_tmp/msout # workaround for CWL bug
+      - id: image_name
+        source: targetname
+        valueFrom: '$(self+"_lowres")'
+      - id: image_size
+        source: getFWHM/fwhm
+        valueFrom: '$((Math.round(2.1 * self * 3600 / 10.)) % 2 > 0 ? [Math.round(2.1 * self * 3600 / 10.) + 1 , Math.round(2.1 * self * 3600 / 10.) + 1] : [Math.round(2.1 * self * 3600 / 10.), Math.round(2.1 * self * 3600 / 10.)])'
+      - id: image_scale
+        default: '30arcsec'
+      - id: weighting
+        default: 'briggs -0.3'
+      - id: no_model_update
+        default: true
+      - id: niter
+        default: 50000
+      - id: minuv-l
+        default: 30
+      - id: parallel-gridding
+        default: 4
+      - id: maxuvw-m
+        default: 6000
+      - id: taper-gaussian
+        default: '200arcsec'
+      - id: baseline_averaging
+        default: getFWHM/fwhm
+        valueFrom: '$(1.87e3 * 60000. * 2. * Math.PI / (24 * 60 * 60 * Math.round(2.1 * self * 3600 / 10.)) < 1 ? false : 1.87e3 * 60000. * 2. * Math.PI / (24 * 60 * 60 * Math.round(2.1 * self * 3600 / 10.)) > 10 ? 10 : 1.87e3 * 60000. * 2. * Math.PI / (24 * 60 * 60 * Math.round(2.1 * self * 3600 / 10.)))'
+      - id: mgain
+        default: 0.85
+      - id: parallel-deconvolution
+        default: 512
+      - id: auto_threshold
+        default: 1.5
+      - id: local_rms
+        default: true
+      - id: auto_mask
+        default: 3.0
+      - id: fits_mask
+        source: blank_image_reg_2/output_image
+      - id: join-channels
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/2.e6) > 1)'
+      - id: channels-out
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/2.e6) > 0 ? Math.round(self/2.e6) : 1)'
+      - id: use-wgridder
+        default: true
+      - id: tempdir
+        source: wsclean_tmpdir
+    out:
+      - id: msout
+      - id: model
+      - id: image
+      - id: logfile
+    run: ../../steps/wsclean.cwl
+    label: image_lowres
+  - id: predict_lowres
+    in:
+      - id: msin
+        source: image_lowres/msout # workaround for CWL bug
+      - id: image_name
+        source: targetname
+        valueFrom: '$(self+"_lowres")'
+      - id: image_size
+        source: getFWHM/fwhm
+        valueFrom: '$((Math.round(2.1 * self * 3600 / 10.)) % 2 > 0 ? [Math.round(2.1 * self * 3600 / 10.) + 1 , Math.round(2.1 * self * 3600 / 10.) + 1] : [Math.round(2.1 * self * 3600 / 10.), Math.round(2.1 * self * 3600 / 10.)])'
+      - id: image_scale
+        default: '30arcsec'
+      - id: fits_image
+        source: image_lowres/model
+      - id: channels-out
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/2.e6) > 0 ? Math.round(self/2.e6) : 1)'
+      - id: use-wgridder
+        default: true
+      - id: tempdir
+        source: wsclean_tmpdir
+    out:
+      - id: msout
+      - id: logfile
+    run: ../../steps/wsclean_predict.cwl
+    label: predict_lowres
+  - id: subtract_lowres
+    in:
+      - id: msin
+        source: predict_lowres/msout
+      - id: command
+        default: "set CORRECTED_DATA = CORRECTED_DATA - MODEL_DATA"
+    out:
+      - id: msout
+      - id: logfile
+    run: ../../steps/taql.cwl
+    label: subtract_lowres
+    scatter:
+      - msin
+  - id: aoflag_residual
+    in:
+      - id: msin
+        source: subtract_lowres/msout
+      - id: verbose
+        default: true
+      - id: concatenate-frequency
+        source: aoflag_freqconcat
+      - id: strategy
+        source: rfistrategy
+      - id: reorder
+        source: aoflag_reorder
+      - id: chunk-size
+        source: aoflag_chunksize
+    out:
+      - id: output_ms
+      - id: logfile
+    run: ../../steps/aoflag.cwl
+    label: aoflag_residual
+  - id: image_large
+    in:
+      - id: msin
+        source: aoflag_residual/output_ms
+      - id: image_name
+        source: targetname
+        valueFrom: '$(self+"_large")'
+      - id: image_size
+        default:
+          - 2000
+          - 2000
+      - id: image_scale
+        default: '20arcsec'
+      - id: no_fit_beam
+        default: true
+      - id: circular_beam
+        default: true
+      - id: beam_size
+        default: 180.0
+      - id: multiscale
+        default: true
+      - id: multiscale_scales
+        default: '0,4,8,16,32,64'
+      - id: weighting
+        default: 'briggs -0.3'
+      - id: no_model_update
+        default: true
+      - id: niter
+        default: 10000
+      - id: minuv-l
+        default: 20
+      - id: maxuvw-m
+        default: 5000
+      - id: taper-gaussian
+        default: '180arcsec'
+      - id: baseline_averaging
+        default: getFWHM/fwhm
+        valueFrom: '$(1.87e3 * 60000. * 2. * Math.PI / (24 * 60 * 60 * Math.round(2.1 * self * 3600 / 10.)) < 1 ? false : 1.87e3 * 60000. * 2. * Math.PI / (24 * 60 * 60 * Math.round(2.1 * self * 3600 / 10.)) > 10 ? 10 : 1.87e3 * 60000. * 2. * Math.PI / (24 * 60 * 60 * Math.round(2.1 * self * 3600 / 10.)))'
+      - id: mgain
+        default: 0.85
+      - id: parallel-deconvolution
+        default: 512
+      - id: auto_threshold
+        default: 0.5
+      - id: local_rms
+        default: true
+      - id: auto_mask
+        default: 1.5
+      - id: join-channels
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/4.e6) > 1)'
+      - id: channels-out
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/4.e6) > 0 ? Math.round(self/4.e6) : 1)'
+      - id: use-wgridder
+        default: true
+      - id: tempdir
+        source: wsclean_tmpdir
+    out:
+      - id: image
+      - id: msout
+      - id: logfile
+    run: ../../steps/wsclean.cwl
+    label: image_large
+  - id: corrupt_model
+    in:
+      - id: max_dp3_threads
+        source: max_dp3_threads
+      - id: msin
+        source: image_large/msout
+      - id: parmdb
+        source: parmdb
+      - id: skymodel_source
+        source: skymodel_source
+    out:
+      - id: msout
+      - id: logfile
+    scatter:
+      - msin
+    run: ./corrupt_model.cwl
+    label: corrupt_model
+  - id: recreate_model
+    in:
+      - id: msin
+        source: corrupt_model/msout
+      - id: image_name
+        source: targetname
+        valueFrom: '$(self+"_mask")'
+      - id: image_size
+        source: getFWHM/fwhm
+        valueFrom: '$((Math.round(2.1 * self * 3600 / 10.)) % 2 > 0 ? [Math.round(2.1 * self * 3600 / 10.) + 1 , Math.round(2.1 * self * 3600 / 10.) + 1] : [Math.round(2.1 * self * 3600 / 10.), Math.round(2.1 * self * 3600 / 10.)])'
+      - id: image_scale
+        default: '10arcsec'
+      - id: fits_model
+        source: image_predict/model
+      - id: channels-out
+        source: total_bandwidth
+        valueFrom: '$(Math.round(self/4.e6) > 0 ? Math.round(self/4.e6) : 1)'
+      - id: use-wgridder
+        default: true
+      - id: tempdir
+        source: wsclean_tmpdir
+    out:
+      - id: msout
+      - id: logfile
+    run: ../../steps/wsclean_predict.cwl
+    label: recreate_model
+  - id: concat_logfiles_predict
+    in:
+      - id: file_list
+        source:
+          - make_mask_image/logfile
+          - make_mask/logfile
+          - blank_image_reg/logfile
+          - image_predict/logfile
+        linkMerge: merge_flattened
+        pickValue: all_non_null
+      - id: file_prefix
+        default: selfcal_predict
+    out:
+      - id: output
+    run: ../../steps/concatenate_files.cwl
+    label: concat_logfiles_makemask
+  - id: merge_array_files
+    in:
+      - id: input
+        source: subtract_model/logfile
+    out:
+      - id: output
+    run: ../../steps/merge_array_files.cwl
+    label: merge_array_files
+  - id: merge_array_files_lowres
+    in:
+      - id: input
+        source: subtract_lowres/logfile
+    out:
+      - id: output
+    run: ../../steps/merge_array_files.cwl
+    label: merge_array_files_lowres
+  - id: concat_logfiles_subtract
+    in:
+      - id: file_list
+        source:
+          - merge_array_files/output
+        linkMerge: merge_flattened
+      - id: file_prefix
+        default: selfcal_subtract
+    out:
+      - id: output
+    run: ../../steps/concatenate_files.cwl
+    label: concat_logfiles_subtract
+  - id: concat_logfiles_lowres
+    in:
+      - id: file_list
+        source:
+          - make_beam_reg/logfile
+          - blank_image_reg_1/logfile
+          - blank_image_reg_2/logfile
+          - image_lowres/logfile
+          - predict_lowres/logfile
+          - merge_array_files_lowres/output
+        linkMerge: merge_flattened
+      - id: file_prefix
+        default: selfcal_lowres
+    out:
+      - id: output
+    run: ../../steps/concatenate_files.cwl
+    label: concat_logfiles_lowres
+  - id: concat_logfiles_aoflag
+    in:
+      - id: file_list
+        source:
+          - aoflag_residual/logfile
+        linkMerge: merge_flattened
+      - id: file_prefix
+        default: aoflag_residual
+    out:
+      - id: output
+    run: ../../steps/concatenate_files.cwl
+    label: concat_logfiles_aoflag
+  - id: concat_logfiles_large
+    in:
+      - id: file_list
+        source:
+          - image_large/logfile
+        linkMerge: merge_flattened
+      - id: file_prefix
+        default: selfcal_large
+    out:
+      - id: output
+    run: ../../steps/concatenate_files.cwl
+    label: concat_logfiles_large
+  - id: concat_logfiles_corrupt
+    in:
+      - id: file_list
+        source:
+          - corrupt_model/logfile
+        linkMerge: merge_flattened
+      - id: file_prefix
+        default: corrupt_model
+    out:
+      - id: output
+    run: ../../steps/concatenate_files.cwl
+    label: concat_logfiles_corrupt
+  - id: concat_logfiles_recreate
+    in:
+      - id: file_list
+        source:
+          - recreate_model/logfile
+        linkMerge: merge_flattened
+      - id: file_prefix
+        default: recreate_model
+    out:
+      - id: output
+    run: ../../steps/concatenate_files.cwl
+    label: concat_logfiles_recreate
+requirements:
+  - class: SubworkflowFeatureRequirement
+  - class: ScatterFeatureRequirement
\ No newline at end of file
diff --git a/workflows/linc_target/prep.cwl b/workflows/linc_target/prep.cwl
index d250f4b1..63cc4724 100644
--- a/workflows/linc_target/prep.cwl
+++ b/workflows/linc_target/prep.cwl
@@ -13,6 +13,12 @@ inputs:
   - id: process_baselines_target
     type: string?
     default: '[CR]S*&'
+  - id: num_SBs_per_group
+    type: int?
+    default: 10
+  - id: reference_stationSB
+    type: int?
+    default: null
   - id: filter_baselines
     type: string?
     default: '[CR]S*&'
@@ -123,6 +129,8 @@ inputs:
   - id: lbfgs_robustdof
     type: float?
     default: 200
+  - id: skymodel_fluxlimit
+    type: float?
 outputs:
   - id: compare_stations_filter
     outputSource:
@@ -151,6 +159,7 @@ outputs:
   - id: logfiles
     outputSource:
       - find_skymodel_target/logfile
+      - sort_times_into_freqGroups/logfile
       - concat_logfiles_skymodels/output
       - check_ateam_separation/logfile
       - check_demix/logfile
@@ -176,7 +185,57 @@ outputs:
     type:
       - File
       - Directory
+  - id: total_bandwidth
+    outputSource:
+      - sort_times_into_freqGroups/total_bandwidth
+    type: float
+  - id: midfreq
+    outputSource:
+      - sort_times_into_freqGroups/midfreq
+    type: float
+  - id: groupnames
+    outputSource:
+      - sort_times_into_freqGroups/groupnames
+    type: string[]
+  - id: filenames
+    outputSource:
+      - sort_times_into_freqGroups/filenames
+    type: File
 steps:
+  - id: sort_times_into_freqGroups
+    in:
+      - id: msin
+        source:
+          - msin
+      - id: numbands
+        source: num_SBs_per_group
+      - id: DP3fill
+        default: true
+      - id: stepname
+        default: .dp3concat
+      - id: firstSB
+        source: reference_stationSB
+      - id: truncateLastSBs
+        default: false
+    out:
+      - id: filenames
+      - id: groupnames
+      - id: total_bandwidth
+      - id: midfreq
+      - id: minfreq
+      - id: logfile
+    run: ../../steps/sort_times_into_freqGroups.cwl
+    label: sorttimesintofreqGroups
+  - id: getFWHM
+    in:
+      - id: msin
+        source: msin
+      - id: beamfreq
+        source: sort_times_into_freqGroups/minfreq
+    out:
+      - id: fwhm
+    run: ../../steps/getFWHM.cwl
+    label: getFWHM
   - id: find_skymodel_target
     in:
       - id: msin
@@ -185,13 +244,16 @@ steps:
       - id: SkymodelPath
         source: target_skymodel
       - id: Radius
-        default: 5
+        source: getFWHM/fwhm
+        valueFrom: $(self/2.0)
       - id: Source
         source: skymodel_source
       - id: DoDownload
         source: use_target
       - id: targetname
         source: targetname
+      - id: fluxlimit
+        source: skymodel_fluxlimit
     out:
       - id: skymodel
       - id: logfile
diff --git a/workflows/linc_target/selfcal_targ.cwl b/workflows/linc_target/selfcal_targ.cwl
index 6b35ef7a..a1e62847 100644
--- a/workflows/linc_target/selfcal_targ.cwl
+++ b/workflows/linc_target/selfcal_targ.cwl
@@ -33,6 +33,33 @@ inputs:
     default: '[CR]S*&'
   - id: insolutions
     type: File
+  - id: targetname
+    type: string?
+    default: 'pointing'
+  - id: wsclean_tmpdir
+    type: string?
+    default: '/tmp'
+  - id: total_bandwidth
+    type: float?
+    default: 50
+  - id: midfreq
+    type: float?
+    default: 50e6
+  - id: selfcal_region
+    type: File?
+  - id: rfistrategy
+    type:
+      - File?
+      - string?
+  - id: aoflag_reorder
+    type: boolean?
+    default: false
+  - id: aoflag_chunksize
+    type: int?
+    default: 2000
+  - id: aoflag_freqconcat
+    type: boolean?
+    default: true
   - id: execute
     type: boolean?
     default: true
@@ -53,12 +80,14 @@ outputs:
     outputSource:
       - fr/logfile
       - tec_and_amp/logfile
+      - imaging_subtract/logfile
     linkMerge: merge_flattened
     type: File[]
   - id: inspection
     outputSource:
       - fr/inspection
       - tec_and_amp/inspection
+      - imaging_subtract/inspection
     type: File[]
     linkMerge: merge_flattened
 steps:
@@ -209,6 +238,40 @@ steps:
     scatter:
       - msin
     label: apply_targ_amp
+  - id: imaging_subtract
+    in:
+      - id: max_dp3_threads
+        source: max_dp3_threads
+      - id: msin
+        source: apply_targ_amp/msout
+      - id: targetname
+        source: targetname
+      - id: wsclean_tmpdir
+        source: wsclean_tmpdir
+      - id: total_bandwidth
+        source: total_bandwidth
+      - id: midfreq
+        source: midfreq
+      - id: selfcal_region
+        source: selfcal_region
+      - id: rfistrategy
+        source: rfistrategy
+      - id: aoflag_reorder
+        source: aoflag_reorder
+      - id: aoflag_chunksize
+        source: aoflag_chunksize
+      - id: aoflag_freqconcat
+        source: aoflag_freqconcat
+      - id: parmdb
+        source: tec_and_amp/outsolutions
+      - id: skymodel_source
+        source: skymodel_source
+    out:
+      - id: msout
+      - id: logfile
+      - id: inspection
+    run: ./imaging_subtract.cwl
+    label: imaging_subtract
 requirements:
   - class: ScatterFeatureRequirement
   - class: SubworkflowFeatureRequirement
\ No newline at end of file
-- 
GitLab