diff --git a/scripts/gather_model_images.sh b/scripts/gather_model_images.sh
new file mode 100755
index 0000000000000000000000000000000000000000..9555938ec0c64d560c76904f06a7d4403423aa06
--- /dev/null
+++ b/scripts/gather_model_images.sh
@@ -0,0 +1,45 @@
+#!/bin/bash
+
+# Ensure a model image folder is provided
+if [ -z "$1" ]; then
+    echo "Usage: $0 <model_image_folder>"
+    exit 1
+fi
+
+MODEL_IMAGE_FOLDER=$1
+
+# Patterns to check with model images
+PATTERNS=(
+    "*-model-fpb.fits"
+    "*-model-pb.fits"
+    "*-model.fits"
+)
+
+# Function to copy files if more than one match exists
+copy_files() {
+    local pattern="$1"
+    local files=()
+
+    # Filter files matching the pattern but excluding those with "MFS"
+    for file in "$MODEL_IMAGE_FOLDER"/$pattern; do
+        [[ "$file" == *"MFS"* ]] && continue  # Skip files with "MFS"
+        files+=("$file")
+    done
+
+    mkdir -p output_images
+
+    if [ "${#files[@]}" -gt 1 ]; then
+        cp "${files[@]}" output_images/
+        return 0
+    fi
+    return 1  # No valid files were copied
+}
+
+# Iterate through patterns and copy the first match
+for pattern in "${PATTERNS[@]}"; do
+    copy_files "$pattern" && exit 0
+done
+
+# If no files were copied, exit with error
+echo "ERROR: missing model images in folder $MODEL_IMAGE_FOLDER"
+exit 1
diff --git a/steps/dp3_parset.cwl b/steps/dp3_parset.cwl
index d44ae6887a8d946ca78371b0a9c61f54985893ba..2672651b561eb758a35e959bff506c7558a510f2 100644
--- a/steps/dp3_parset.cwl
+++ b/steps/dp3_parset.cwl
@@ -2,8 +2,7 @@ cwlVersion: v1.2
 class: CommandLineTool
 id: dp3_parset
 label: DP3 with parset
-doc: |
-    Run DP3 with a parset
+doc: Run DP3 with a parset
 
 baseCommand: DP3
 
@@ -22,7 +21,7 @@ outputs:
     type: Directory
     doc: Output measurement set
     outputBinding:
-      glob: "*.ms"
+      glob: "*.concat.ms"
   - id: logfile
     type: File[]
     outputBinding:
diff --git a/steps/gather_model_images.cwl b/steps/gather_model_images.cwl
new file mode 100644
index 0000000000000000000000000000000000000000..669e89c871c5398147c1f12dd7751a82dd2b3b8e
--- /dev/null
+++ b/steps/gather_model_images.cwl
@@ -0,0 +1,25 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: gatherdis2
+doc: Gather the correct WSClean model images from model_image_folder
+
+baseCommand:
+  - gather_model_images.sh
+
+inputs:
+  - id: model_image_folder
+    type: Directory
+    doc: Directory with model images
+    inputBinding:
+      position: 1
+
+outputs:
+  - id: filtered_model_image_folder
+    type: Directory
+    doc: Directory with filtered WSClean model images
+    outputBinding:
+      glob: output_images
+
+hints:
+  - class: DockerRequirement
+    dockerPull: vlbi-cwl
diff --git a/steps/get_facet_layout.cwl b/steps/get_facet_layout.cwl
new file mode 100644
index 0000000000000000000000000000000000000000..aabbbba11614c35d513f520bbb3687447184fff5
--- /dev/null
+++ b/steps/get_facet_layout.cwl
@@ -0,0 +1,81 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: get_facet_layout
+label: DS9 Facet Generator
+doc: Generates DS9 facet layout for direction-dependent facet imaging.
+
+baseCommand: python3
+
+inputs:
+  - id: msin
+    type: Directory
+    doc: MeasurementSet
+    inputBinding:
+      prefix: "--ms"
+      position: 2
+
+  - id: h5parm
+    type: File
+    doc: Multi-directional HDF5 file.
+    inputBinding:
+      prefix: "--h5"
+      position: 3
+
+  - id: imsize
+    type: int
+    doc: Image size in pixels (larger than image fov to help wsclean 1.2" imaging at boundaries).
+    default: 25000
+    inputBinding:
+      prefix: "--imsize"
+      position: 4
+
+  - id: pixelscale
+    type: float
+    doc: Pixel scale in arcseconds per pixel.
+    default: 0.4
+    inputBinding:
+      prefix: "--pixelscale"
+      position: 5
+
+  - id: regionfile
+    type: string
+    doc: Name of the output DS9 region file.
+    default: "facets.reg"
+    inputBinding:
+      prefix: "--DS9regionout"
+      position: 6
+
+  - id: facetselfcal
+    type: Directory
+    doc: facetselfcal directory.
+
+outputs:
+  - id: facet_regions
+    type: File
+    doc: The output DS9 region file.
+    outputBinding:
+      glob: $(inputs.regionfile)
+  - id: logfile
+    type: File[]
+    doc: log files from get_facet_layout
+    outputBinding:
+      glob: get_facet_layout*.log
+
+arguments:
+  - $( inputs.facetselfcal.path + '/submods/ds9facetgenerator.py' )
+
+requirements:
+  - class: InlineJavascriptRequirement
+  - class: InitialWorkDirRequirement
+    listing:
+      - entry: $(inputs.h5parm)
+      - entry: $(inputs.msin)
+
+hints:
+  - class: DockerRequirement
+    dockerPull: vlbi-cwl
+  - class: ResourceRequirement
+    coresMin: 1
+
+stdout: get_facet_layout.log
+stderr: get_facet_layout_err.log
diff --git a/steps/make_concat_parsets.cwl b/steps/make_concat_parsets.cwl
index db297165d13d45223784aefd19e18166779d2139..f1586f1b8113f5fdcff4cdc24461b0701c0a55a5 100644
--- a/steps/make_concat_parsets.cwl
+++ b/steps/make_concat_parsets.cwl
@@ -16,6 +16,14 @@ inputs:
         position: 1
         separate: true
     doc: Input data in MeasurementSet format.
+  - id: dysco_bitrate
+    type: int?
+    doc: Number of bits per float used for columns containing visibilities.
+    default: 8
+    inputBinding:
+        prefix: "--bitrate"
+        position: 2
+        separate: true
   - id: lofar_helpers
     type: Directory
     doc: Path to lofar_helpers directory.
diff --git a/steps/predict_facet.cwl b/steps/predict_facet.cwl
new file mode 100644
index 0000000000000000000000000000000000000000..ec54628581cfa1f1bf1f31955a3b0c70452ff68d
--- /dev/null
+++ b/steps/predict_facet.cwl
@@ -0,0 +1,110 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: predict_facet
+label: Predict with WSClean
+doc: Uses WSClean to predict sources within a facet and adds the predicted visibilities to the input data.
+
+
+baseCommand: python3
+
+inputs:
+  - id: subtracted_ms
+    type: Directory
+    doc: Input data in MeasurementSet format.
+    inputBinding:
+      prefix: "--mslist"
+      position: 1
+
+  - id: model_image_folder
+    type: Directory
+    doc: Folder containing 1.2" model images.
+    inputBinding:
+      prefix: "--model_image_folder"
+      position: 2
+
+  - id: polygon_region
+    type: File
+    doc: DS9 region file with facets for prediction.
+    inputBinding:
+      prefix: "--region"
+      position: 3
+
+  - id: h5parm
+    type: File
+    doc: HDF5 file containing the solutions for prediction.
+    inputBinding:
+      prefix: "--h5parm_predict"
+      position: 4
+
+  - id: lofar_helpers
+    type: Directory
+    doc: LOFAR helpers directory.
+
+  - id: polygon_info
+    type: File
+    doc: CSV file with polygon information (RA/DEC of calibrator and facet centers and averaging factor)
+
+  - id: copy_to_local_scratch
+    type: boolean?
+    default: false
+    inputBinding:
+      prefix: "--copy_to_local_scratch"
+      position: 5
+      separate: false
+    doc: Specific option for using --bypass-file-store on the Spider cluster to run predict and subtract on local scratch.
+
+  - id: ncpu
+    type: int?
+    doc: Number of cores to use during the subtract.
+    default: 16
+
+outputs:
+  - id: logfile
+    type: File[]
+    doc: Log files from current step.
+    outputBinding:
+      glob: predict_facet*.log
+  - id: facet_ms
+    type: Directory
+    doc: MeasurementSet after predicting back facet model visibilities
+    outputBinding:
+      glob: facet*.ms
+
+
+arguments:
+  - $( inputs.lofar_helpers.path + '/subtract/subtract_with_wsclean.py' )
+  - --applybeam
+  - --applycal
+  - --forwidefield
+  - --inverse
+  - --speedup_facet_subtract
+  - --cleanup_input_ms
+
+requirements:
+  - class: InlineJavascriptRequirement
+  - class: InitialWorkDirRequirement
+    listing: >
+      ${
+        // Set 'writable' on the "subtracted_ms" entry only if copy_to_local_scratch is true.
+        let stagedListing = [
+          { entry: inputs.subtracted_ms },
+          { entry: inputs.model_image_folder },
+          { entry: inputs.polygon_region },
+          { entry: inputs.h5parm },
+          { entry: inputs.polygon_info }
+        ];
+        if (!inputs.copy_to_local_scratch) {
+          stagedListing[0].writable = true;
+        }
+        return stagedListing;
+      }
+
+
+hints:
+  - class: DockerRequirement
+    dockerPull: vlbi-cwl
+  - class: ResourceRequirement
+    coresMin: $(inputs.ncpu)
+
+stdout: predict_facet.log
+stderr: predict_facet_err.log
\ No newline at end of file
diff --git a/steps/remove_flagged_stations.cwl b/steps/remove_flagged_stations.cwl
new file mode 100644
index 0000000000000000000000000000000000000000..68b23fda4c00053576dd1b2cd822dd4c70d4b178
--- /dev/null
+++ b/steps/remove_flagged_stations.cwl
@@ -0,0 +1,51 @@
+cwlVersion: v1.2
+class: CommandLineTool
+id: remove_flagged_stations
+label: Remove fully flagged stations
+doc: Removes from the MeasurementSet all stations for which all data are flagged.
+
+baseCommand:
+  - python3
+
+inputs:
+    - id: ms
+      type: Directory
+      doc: MeasurementSet
+      inputBinding:
+        position: 3
+    - id: lofar_helpers
+      type: Directory
+      doc: LOFAR helpers directory.
+
+outputs:
+    - id: cleaned_ms
+      type: Directory
+      doc: MeasurementSet where fully flagged stations are removed.
+      outputBinding:
+        glob: $( 'flagged_' + inputs.ms.basename )
+    - id: logfile
+      type: File[]
+      doc: Log files from current step.
+      outputBinding:
+        glob: remove_flagged_stations*.log
+
+requirements:
+  - class: InlineJavascriptRequirement
+  - class: InitialWorkDirRequirement
+    listing:
+      - entry: $(inputs.ms)
+        writable: false
+
+arguments:
+  - $( inputs.lofar_helpers.path + '/ms_helpers/remove_flagged_stations.py' )
+  - --msout
+  - $( 'flagged_' + inputs.ms.basename )
+
+hints:
+  - class: DockerRequirement
+    dockerPull: vlbi-cwl
+  - class: ResourceRequirement
+    coresMin: 2
+
+stdout: remove_flagged_stations.log
+stderr: remove_flagged_stations_err.log
\ No newline at end of file
diff --git a/steps/split_polygons.cwl b/steps/split_polygons.cwl
new file mode 100644
index 0000000000000000000000000000000000000000..9b71e4e0f0fdd0a3bd6f60214ac049f7c92c9d98
--- /dev/null
+++ b/steps/split_polygons.cwl
@@ -0,0 +1,63 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: split_polygons
+label: Split Polygon Facets
+doc: This step splits a multi-facet region file into individual facet region files.
+
+baseCommand: python3
+
+inputs:
+  - id: h5parm
+    type: File
+    doc: Multi-directional HDF5 file.
+    inputBinding:
+      prefix: "--h5"
+      position: 1
+
+  - id: facet_regions
+    type: File
+    doc: The DS9 region file that defines the facets.
+    inputBinding:
+      prefix: "--reg"
+      position: 2
+
+  - id: lofar_helpers
+    type: Directory
+    doc: The lofar_helpers directory.
+
+outputs:
+  - id: polygon_regions
+    type: File[]
+    doc: The facet regions.
+    outputBinding:
+      glob: "poly*.reg"
+  - id: polygon_info
+    type: File
+    doc: Polygon CSV file.
+    outputBinding:
+      glob: "*.csv"
+  - id: logfile
+    type: File[]
+    doc: Log files from current step.
+    outputBinding:
+      glob: split_polygons*.log
+
+arguments:
+  - $(inputs.lofar_helpers.path)/ds9_helpers/split_polygon_facets.py
+  - --extra_boundary=0.0
+
+requirements:
+  - class: InlineJavascriptRequirement
+  - class: InitialWorkDirRequirement
+    listing:
+      - entry: $(inputs.h5parm)
+      - entry: $(inputs.facet_regions)
+
+hints:
+  - class: DockerRequirement
+    dockerPull: vlbi-cwl
+  - class: ResourceRequirement
+    coresMin: 1
+
+stdout: split_polygons.log
+stderr: split_polygons_err.log
diff --git a/steps/subtract_fov_wsclean.cwl b/steps/subtract_fov_wsclean.cwl
new file mode 100644
index 0000000000000000000000000000000000000000..3df55740b4d5b9569f4e8b8f74f9fd3daee25a81
--- /dev/null
+++ b/steps/subtract_fov_wsclean.cwl
@@ -0,0 +1,97 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: subtract_with_wsclean
+label: Subtract with WSClean
+doc: This step uses WSClean to subtract visibilities corresponding to model images.
+
+baseCommand: python3
+
+inputs:
+  - id: msin
+    type: Directory
+    doc: MeasurementSet for source subtraction.
+    inputBinding:
+      prefix: "--mslist"
+      position: 1
+
+  - id: model_image_folder
+    type: Directory
+    doc: Directory containing 1.2" (or optionally other resolution) model images.
+    inputBinding:
+      prefix: "--model_image_folder"
+      position: 2
+
+  - id: facet_regions
+    type: File
+    doc: The DS9 region file that defines the facets for prediction.
+    inputBinding:
+      prefix: "--facets_predict"
+      position: 3
+
+  - id: h5parm
+    type: File
+    doc: The HDF5 solution file containing the solutions for prediction.
+    inputBinding:
+      prefix: "--h5parm_predict"
+      position: 4
+
+  - id: lofar_helpers
+    type: Directory
+    doc: LOFAR helpers directory.
+
+  - id: copy_to_local_scratch
+    type: boolean?
+    default: false
+    inputBinding:
+      prefix: "--copy_to_local_scratch"
+      position: 5
+      separate: false
+    doc: Specific option for using --bypass-file-store on the Spider cluster to run predict and subtract on local scratch.
+
+  - id: ncpu
+    type: int?
+    doc: Number of cores to use during the subtract.
+    default: 16
+
+outputs:
+  - id: logfile
+    type: File[]
+    doc: Log files from current step.
+    outputBinding:
+      glob: subtract_fov*.log
+  - id: subtracted_ms
+    type: Directory
+    doc: MS subtracted data
+    outputBinding:
+      glob: subfov_*.ms
+
+
+arguments:
+  - $( inputs.lofar_helpers.path + '/subtract/subtract_with_wsclean.py' )
+
+requirements:
+  - class: InlineJavascriptRequirement
+  - class: InitialWorkDirRequirement
+    listing: >
+      ${
+        // Set 'writable' on the "msin" entry only if copy_to_local_scratch is true.
+        let stagedListing = [
+          { entry: inputs.msin },
+          { entry: inputs.model_image_folder },
+          { entry: inputs.facet_regions },
+          { entry: inputs.h5parm }
+        ];
+        if (!inputs.copy_to_local_scratch) {
+          stagedListing[0].writable = true;
+        }
+        return stagedListing;
+      }
+
+hints:
+  - class: DockerRequirement
+    dockerPull: vlbi-cwl
+  - class: ResourceRequirement
+    coresMin: $(inputs.ncpu)
+
+stdout: subtract_fov.log
+stderr: subtract_fov_err.log
\ No newline at end of file
diff --git a/workflows/facet_subtract.cwl b/workflows/facet_subtract.cwl
new file mode 100644
index 0000000000000000000000000000000000000000..7bd83f550b3ca794cf98aff87e7f3277dece83ba
--- /dev/null
+++ b/workflows/facet_subtract.cwl
@@ -0,0 +1,162 @@
+class: Workflow
+cwlVersion: v1.2
+id: facet_subtract
+label: Facet subtraction
+doc: Use WSClean to predict and subtract model data, to split all facets into separate MeasurementSets.
+
+inputs:
+    - id: msin
+      type: Directory[]
+      doc: Unaveraged MeasurementSets with coverage of the target directions.
+    - id: model_image_folder
+      type: Directory
+      doc: Folder with 1.2" model images.
+    - id: h5parm
+      type: File
+      doc: Merged h5parms
+    - id: lofar_helpers
+      type: Directory
+      doc: LOFAR helpers directory.
+    - id: facetselfcal
+      type: Directory
+      doc: facetselfcal directory.
+    - id: copy_to_local_scratch
+      type: boolean?
+      doc: Specific option for using --bypass-file-store on the Spider cluster to run predict and subtract on local scratch.
+      default: false
+    - id: ncpu
+      type: int?
+      doc: Number of cores to use during predict and subtract.
+      default: 16
+    - id: dysco_bitrate
+      type: int?
+      doc: Number of bits per float used for columns containing visibilities.
+      default: 8
+
+steps:
+    - id: get_facet_layout
+      label: Get DS9 facet layout
+      in:
+        - id: msin
+          source: msin
+          valueFrom: $(self[0])
+        - id: h5parm
+          source: h5parm
+        - id: facetselfcal
+          source: facetselfcal
+      out:
+        - id: facet_regions
+      run: ../steps/get_facet_layout.cwl
+
+    - id: get_model_images
+      label: Get WSClean model images
+      in:
+        - id: model_image_folder
+          source: model_image_folder
+      out:
+        - id: filtered_model_image_folder
+      run: ../steps/gather_model_images.cwl
+
+    - id: subtract_fov_wsclean
+      label: Subtract complete FoV
+      in:
+         - id: msin
+           source: msin
+         - id: h5parm
+           source: h5parm
+         - id: facet_regions
+           source: get_facet_layout/facet_regions
+         - id: model_image_folder
+           source: get_model_images/filtered_model_image_folder
+         - id: lofar_helpers
+           source: lofar_helpers
+         - id: copy_to_local_scratch
+           source: copy_to_local_scratch
+         - id: ncpu
+           source: ncpu
+      out:
+         - subtracted_ms
+      run: ../steps/subtract_fov_wsclean.cwl
+      scatter: msin
+
+    - id: split_polygons
+      label: Split polygon file
+      in:
+         - id: facet_regions
+           source: get_facet_layout/facet_regions
+         - id: h5parm
+           source: h5parm
+         - id: lofar_helpers
+           source: lofar_helpers
+      out:
+         - id: polygon_info
+         - id: polygon_regions
+      run: ../steps/split_polygons.cwl
+
+    - id: predict_facet
+      label: Predict a polygon back in empty MS
+      in:
+         - id: subtracted_ms
+           source: subtract_fov_wsclean/subtracted_ms
+         - id: polygon_region
+           source: split_polygons/polygon_regions
+         - id: h5parm
+           source: h5parm
+         - id: polygon_info
+           source: split_polygons/polygon_info
+         - id: model_image_folder
+           source: get_model_images/filtered_model_image_folder
+         - id: lofar_helpers
+           source: lofar_helpers
+         - id: copy_to_local_scratch
+           source: copy_to_local_scratch
+         - id: ncpu
+           source: ncpu
+      out:
+         - facet_ms
+      run: ../steps/predict_facet.cwl
+      scatter: [subtracted_ms, polygon_region]
+      scatterMethod: flat_crossproduct
+
+    - id: make_concat_parset
+      label: Make concat parsets
+      in:
+         - id: msin
+           source: predict_facet/facet_ms
+         - id: lofar_helpers
+           source: lofar_helpers
+         - id: dysco_bitrate
+           source: dysco_bitrate
+      out:
+         - id: concat_parsets
+      run: ../steps/make_concat_parsets.cwl
+
+    - id: concat_facets
+      label: Run DP3 parsets
+      in:
+        - id: parset
+          source: make_concat_parset/concat_parsets
+        - id: msin
+          source: predict_facet/facet_ms
+      out:
+        - id: msout
+      run: ../steps/dp3_parset.cwl
+      scatter: parset
+
+requirements:
+  - class: MultipleInputFeatureRequirement
+  - class: ScatterFeatureRequirement
+
+outputs:
+    - id: facet_ms
+      type: Directory[]
+      outputSource:
+        - concat_facets/msout
+        - predict_facet/facet_ms
+      pickValue: first_non_null
+    - id: polygon_info
+      type: File
+      outputSource: split_polygons/polygon_info
+    - id: polygon_regions
+      type: File[]
+      outputSource: split_polygons/polygon_regions