diff --git a/scripts/make_pol_cubes.py b/scripts/make_pol_cubes.py
new file mode 100755
index 00000000..a0b8e85c
--- /dev/null
+++ b/scripts/make_pol_cubes.py
@@ -0,0 +1,142 @@
+#!/usr/bin/env python3
+"""Create Stokes Q and U cubes from WSClean channel images"""
+
+
+import numpy as np
+from astropy.io import fits
+from typing import List
+import os
+import argparse
+
+from validate_lofar_images import get_rms
+
+
+def cube_maker(q_images: List[str], u_images: List[str], nchan: int, imsize: List[int]) -> None:
+    """
+    Constructs Stokes Q and U data cubes from individual FITS channel images from WSClean,
+    filtering out noisy channels based on an RMS threshold.
+    This function reads a sequence of Stokes Q and U FITS files, extracts their
+    data and frequencies, and calculates the average RMS noise per channel. Channels
+    with an average noise greater than 5 times the median noise across all channels
+    are flagged and excluded. The filtered results are written out as new 4D FITS 
+    cubes alongside metadata text files.
+
+    Parameters
+    ----------
+    q_images : list of str
+        File paths to the Stokes Q FITS images, ordered by channel.
+    u_images : list of str
+        File paths to the Stokes U FITS images, ordered by channel.
+    nchan : int
+        The total number of frequency channels to process.
+    imsize : tuple of int or list of int
+        The spatial dimensions of the input images in pixels.
+
+    Returns
+    -------
+    None
+        Outputs are written directly to disk as FITS files 
+        ('<imagename>-polcube-Q.fits', '<imagename>-polcube-U.fits')
+        and text data files.
+
+    Notes
+    -----
+    - The output files use a naming prefix extracted from the first Stokes Q filename.
+    - Temporary logs of frequencies and noise values are generated and then 
+        overwritten with the final filtered subsets.
+    """
+
+    # Extract the name for the frequency and rms noise files
+    firstq = os.path.basename(q_images[0])
+    imagename = firstq.split('-')[0]
+
+    cube_q = np.zeros((1,nchan,imsize[1],imsize[0]))
+    cube_u = np.zeros((1,nchan,imsize[1],imsize[0]))
+
+    header_q = None
+    header_u = None
+    
+    with open('' + str(imagename) + '_frequency_list.dat','w') as lfr:
+        with open('' + str(imagename) + '_avg_qunoise_list.dat','w') as avgqunoise:
+            for i in range(0,nchan):
+                hdu_q = fits.open(q_images[i])
+                hdu_u = fits.open(u_images[i])
+                data_q = hdu_q[0].data[:,:]
+                data_u = hdu_u[0].data[:,:]
+                if np.isnan(data_q).any() or np.isnan(data_u).any():
+                    print("Channel " + str(i) + "is NaN, move to the next one...")
+                    continue
+
+                if header_q is None:      # This is required to store the first frequency channel in the header of the final cube, useful for visualization
+                    header_q = hdu_q[0].header
+                    header_u = hdu_u[0].header
+
+                avgnoise = 0.5 * ( get_rms(q_images[i]) + get_rms(u_images[i]) )
+                frequ = hdu_q[0].header['CRVAL3']
+                lfr.write(str(frequ)+'\n')
+                avgqunoise.write(str(avgnoise)+'\n')
+                cube_q[0,i,:,:] = data_q
+                cube_u[0,i,:,:] = data_u
+
+
+    # Check the average QU noise for each channel. If it is more than 5 times the
+    # median noise, then it is required to exclude that channel.
+
+    with open('' + str(imagename) + '_avg_qunoise_list.dat','r') as avgqunoise:
+        avgqunoise_values = np.array([float(line.strip()) for line in avgqunoise])
+
+    with open('' + str(imagename) + '_frequency_list.dat','r') as lfr:
+        frequencies = np.array([float(line.strip()) for line in lfr])
+
+    mask = avgqunoise_values <= 5. * np.median(avgqunoise_values)
+
+    correct_frequencies = frequencies[mask]
+    correct_avgqunoise = avgqunoise_values[mask]
+    correct_cube_q = cube_q[:,mask,:,:]
+    correct_cube_u = cube_u[:,mask,:,:]
+
+    with open('' + str(imagename) + '_frequency_list.dat','w') as lfr:
+        for freq in correct_frequencies:
+            lfr.write(str(freq)+'\n')
+    
+    with open('' + str(imagename) + '_avg_qunoise_list.dat','w') as avgqunoise:
+        for noise in correct_avgqunoise:
+            avgqunoise.write(str(noise)+'\n')
+
+    
+    # Writing the cubes
+    hdu_cube_q = fits.PrimaryHDU(correct_cube_q,header_q)
+    hdu_cube_q.writeto('' + str(imagename) + '-polcube-Q.fits', overwrite=True)
+    print("Stokes Q cube written")
+
+    hdu_cube_u = fits.PrimaryHDU(correct_cube_u,header_u)
+    hdu_cube_u.writeto('' + str(imagename) + '-polcube-U.fits', overwrite=True)
+    print("Stokes U cube written")
+
+
+def main():
+    parser = argparse.ArgumentParser(description='Cubes maker from WSClean images for RM-synthesis with RMtools')
+    parser.add_argument('--qimages', help='List of input Stokes Q channel images', type=str, required=True)
+    parser.add_argument('--uimages', help='List of input Stokes U channel images', type=str, required=True)
+    parser.add_argument('--nchan', help='Channels out as in WSClean', default=480, type=int)
+    parser.add_argument('--imsize', help='Image size in pixels as in WSClean (width,height or single value)', default='1024,1024', type=str)
+    args = parser.parse_args()
+
+    q_images = args.qimages.split(',')
+    u_images = args.uimages.split(',')
+
+    q_images = sorted(q_images)
+    u_images = sorted(u_images)
+
+    imsize_parts = args.imsize.split(',')
+    if len(imsize_parts) == 2:
+        imsize = [int(imsize_parts[0]), int(imsize_parts[1])]
+    else:
+        single_size = int(imsize_parts[0])
+        imsize = [single_size, single_size]
+
+    cube_maker(q_images, u_images, args.nchan, imsize)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/scripts/run_rmtools.py b/scripts/run_rmtools.py
new file mode 100755
index 00000000..d0ac5e4d
--- /dev/null
+++ b/scripts/run_rmtools.py
@@ -0,0 +1,96 @@
+#!/usr/bin/env python3
+"""Run RM-Tools rmsynth3d locally."""
+
+import argparse
+import glob
+import os
+import shutil
+import subprocess
+import sys
+from typing import List
+from pathlib import Path
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(
+        description="Run RM-Tools rmsynth3d on Stokes Q/U cubes and a frequency list.",
+    )
+    parser.add_argument("--stokes-q", required=True, help="Stokes Q cube (per-channel).")
+    parser.add_argument("--stokes-u", required=True, help="Stokes U cube (per-channel).")
+    parser.add_argument("--freqs", required=True, help="Frequency list in Hz (one per channel).")
+    parser.add_argument("--l", dest="max_lam2", default="150", help="Maximum lambda-squared (-l).")
+    parser.add_argument("--d", dest="dlam2", default="0.3", help="Lambda-squared channel width (-d).")
+    parser.add_argument(
+        "--output-prefix",
+        dest="output_prefix",
+        default="",
+        help="Prefix to prepend to output files (passed to rmsynth3d -o).",
+    )
+    parser.add_argument(
+        "--extra-args",
+        default="",
+        help="Extra arguments passed to rmsynth3d (as a single string).",
+    )
+    return parser.parse_args()
+
+def stage_inputs(args: argparse.Namespace):
+    workdir = Path.cwd()
+
+    q_local = workdir / "image-polcube-Q.fits"
+    u_local = workdir / "image-polcube-U.fits"
+    freq_local = workdir / "image_frequency_list.dat"
+
+    shutil.copy2(args.stokes_q, q_local)
+    shutil.copy2(args.stokes_u, u_local)
+    shutil.copy2(args.freqs, freq_local)
+
+    return str(q_local), str(u_local), str(freq_local)
+
+def build_rmsynth_cmd(args: argparse.Namespace, qfile: str, ufile: str, freqfile: str,) -> List[str]:
+    cmd = [
+        "rmsynth3d",
+        qfile,
+        ufile,
+        freqfile,
+        "-l",
+        str(args.max_lam2),
+        "-d",
+        str(args.dlam2),
+        "-v",
+        "-R",
+    ]
+    if args.output_prefix:
+        cmd += ["-o", args.output_prefix]
+    if args.extra_args:
+        cmd += args.extra_args.split()
+    return cmd
+
+
+def move_outputs(stokes_q_path: str) -> None:
+    # RM-Tools writes outputs into the same directory as the staged input FITS.
+    input_dir = os.path.dirname(stokes_q_path)
+    moved = []
+
+    for path in glob.glob(os.path.join(input_dir, "*FDF_*.fits")):
+        dest = os.path.basename(path)
+        shutil.move(path, dest)
+        moved.append(dest)
+
+    if moved:
+        print("Moved outputs to workdir:")
+        for name in moved:
+            print(f"  {name}")
+    else:
+        print("No staged outputs found to move.")
+
+
+def main() -> int:
+    args = parse_args()
+    qfile,ufile,freqfile = stage_inputs(args)
+    rmsynth_cmd = build_rmsynth_cmd(args,qfile,ufile,freqfile,)
+    print("Running:", " ".join(rmsynth_cmd))
+    print("Working directory: ", Path.cwd())
+    return subprocess.call(rmsynth_cmd)
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/scripts/validate_1arcsec_image.py b/scripts/validate_1arcsec_image.py
old mode 100644
new mode 100755
diff --git a/steps/concat_pol.cwl b/steps/concat_pol.cwl
new file mode 100644
index 00000000..0633c2c4
--- /dev/null
+++ b/steps/concat_pol.cwl
@@ -0,0 +1,62 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: make_cubes
+doc: Concatenate Q and U images into cubes for RM-synthesis
+
+baseCommand:
+  - make_pol_cubes.py
+
+inputs:
+  - id: Q_images
+    type: File[]
+    doc: List of Stokes Q images from /steps/wsclean_pol.cwl.
+    inputBinding:
+      position: 0
+      prefix: '--qimages'
+      itemSeparator: ','
+  - id: U_images
+    type: File[]
+    doc: List of Stokes U images from /steps/wsclean_pol.cwl.
+    inputBinding:
+      position: 0
+      prefix: '--uimages'
+      itemSeparator: ','
+  - id: image_size
+    type: int[]
+    doc: Image size in pixels as [x, y], matching WSClean -size.
+    inputBinding:
+      position: 0
+      prefix: '--imsize'
+      itemSeparator: ','
+  - id: nchannels
+    type: int
+    doc: Number of channels as provided in WSClean as -channels-out.
+    default: 480
+    inputBinding:
+      position: 0
+      prefix: '--nchan'
+
+outputs:
+  - id: stokesQcube
+    type: File
+    doc: Name of the Stokes Q that will be created.
+    outputBinding:
+      glob: "*-polcube-Q.fits"
+  - id: stokesUcube
+    type: File
+    doc: Name of the Stokes U that will be created.
+    outputBinding:
+      glob: "*-polcube-U.fits"
+  - id: frequencies_list
+    type: File
+    doc: List of channels frequencies.
+    outputBinding:
+      glob: "*frequency_list.dat"
+  - id: rms_list
+    type: File
+    doc: List of average Q and U rms noise per channel.
+    outputBinding:
+      glob: "*avg_qunoise_list.dat"
+
+requirements:
+    - class: InlineJavascriptRequirement
diff --git a/steps/run_rmtools.cwl b/steps/run_rmtools.cwl
new file mode 100644
index 00000000..a8e6fed5
--- /dev/null
+++ b/steps/run_rmtools.cwl
@@ -0,0 +1,93 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: run_rmtools
+label: Run RM-Tools rmsynth3d
+doc: |
+  Run RM-Tools rmsynth3d on Stokes Q/U cubes and a frequency list.
+
+baseCommand: run_rmtools.py
+inputs:
+  - id: stokes_q
+    type: File
+    doc: Stokes Q cube (per-channel).
+    inputBinding:
+      position: 1
+      prefix: --stokes-q
+  - id: stokes_u
+    type: File
+    doc: Stokes U cube (per-channel).
+    inputBinding:
+      position: 1
+      prefix: --stokes-u
+  - id: freqs_hz
+    type: File
+    doc: Frequency list in Hz (one per channel).
+    inputBinding:
+      position: 1
+      prefix: --freqs
+  - id: max_lam2
+    type: float?
+    default: 150
+    doc: Maximum lambda-squared value for rmsynth3d (-l).
+    inputBinding:
+      position: 1
+      prefix: --l
+  - id: dlam2
+    type: float?
+    default: 0.3
+    doc: Lambda-squared channel width for rmsynth3d (-d).
+    inputBinding:
+      position: 1
+      prefix: --d
+  - id: output_prefix
+    type: string?
+    doc: Prefix for output products passed to rmsynth3d (-o). Defaults to Stokes Q basename if not provided.
+    inputBinding:
+      position: 1
+      prefix: --output-prefix
+  - id: extra_args
+    type: string?
+    doc: Extra arguments passed to rmsynth3d.
+    inputBinding:
+      position: 1
+      prefix: --extra-args
+outputs:
+  - id: fdf_im_dirty
+    type: File
+    outputBinding:
+      glob: '*FDF_im_dirty.fits'
+  - id: fdf_real_dirty
+    type: File
+    outputBinding:
+      glob: '*FDF_real_dirty.fits'
+  - id: fdf_tot_dirty
+    type: File
+    outputBinding:
+      glob: '*FDF_tot_dirty.fits'
+  - id: fdf_maxpi
+    type: File
+    outputBinding:
+      glob: '*FDF_maxPI.fits'
+  - id: fdf_peakrm
+    type: File
+    outputBinding:
+      glob: '*FDF_peakRM.fits'
+  - id: rmsynth_stdout
+    type: File
+    outputBinding:
+      glob: rmsynth3d_stdout.log
+  - id: rmsynth_stderr
+    type: File
+    outputBinding:
+      glob: rmsynth3d_stderr.log
+
+requirements:
+  - class: InlineJavascriptRequirement
+hints:
+  - class: DockerRequirement
+    dockerPull: vlbi-cwl
+
+
+stdout: rmsynth3d_stdout.log
+stderr: rmsynth3d_stderr.log
+
diff --git a/steps/wsclean_pol.cwl b/steps/wsclean_pol.cwl
new file mode 100644
index 00000000..1b7d91a7
--- /dev/null
+++ b/steps/wsclean_pol.cwl
@@ -0,0 +1,294 @@
+class: CommandLineTool
+cwlVersion: v1.2
+id: wsclean_pol
+label: WSClean
+doc: Runs WSClean on the input data to produce an image.
+
+baseCommand: wsclean
+arguments: [-verbose, -log-time, -no-update-model-required]
+
+inputs:
+  - id: msin
+    type:
+      - Directory
+      - Directory[]
+    inputBinding:
+      position: 2
+      shellQuote: false
+      itemSeparator: ' '
+  - id: tempdir
+    type: string
+    default: '.'
+    inputBinding:
+      position: 1
+      shellQuote: false
+      itemSeparator: ' '
+      prefix: '-temp-dir'
+  - id: cores
+    type: int?
+    default: 24
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-j'
+  - id: size
+    type: int[]?
+    default: [1024, 1024]
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-size'
+  - id: baseline_averaging
+    type: float?
+    default: 7.74
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-baseline-averaging'
+  - id: minuv-l
+    type: float?
+    default: 80.0
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-minuv-l'
+  - id: weight
+    type:
+      - string?
+    default: briggs -1.5
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-weight'
+  - id: parallel-reordering
+    type: int?
+    default: 6
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-parallel-reordering'
+  - id: mgain
+    type: float?
+    default: 0.7
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-mgain'
+  - id: data-column
+    type: string?
+    default: CORRECTED_DATA
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-data-column'
+  - id: auto-mask
+    type: float?
+    default: 3.0
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-auto-mask'
+  - id: auto-threshold
+    type: float?
+    default: 1.0
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-auto-threshold'
+  - id: pol
+    type: string?
+    default: i,q,u,v
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-pol'
+  - id: name
+    type: string?
+    default: "image"
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-name'
+  - id: scale
+    type: string?
+    default: "0.75asec"
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-scale'
+  - id: beam-size
+    type: string?
+    default: "0.3asec"
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-beam-size'
+  - id: taper-gaussian
+    type: string?
+    default: 1.2asec
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-taper-gaussian'
+  - id: niter
+    type: int
+    default: 150000
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-niter'
+  - id: multiscale-scale-bias
+    type: float?
+    default: 0.6
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-multiscale-scale-bias'
+  - id: parallel-deconvolution
+    type: int?
+    default: 2600
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-parallel-deconvolution'
+  - id: parallel-gridding
+    type: int?
+    default: 4
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-parallel-gridding'
+  - id: multiscale
+    type: boolean?
+    default: true
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-multiscale'
+  - id: multiscale-max-scales
+    type: int?
+    default: 3
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-multiscale-max-scales'
+  - id: nmiter
+    type: int?
+    default: 9
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-nmiter'
+  - id: channels-out
+    type: int?
+    default: 480
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-channels-out'
+  - id: join-channels
+    type: boolean?
+    default: true
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-join-channels'
+  - id: join-polarizations
+    type: boolean?
+    default: true
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-join-polarizations'
+  - id: squared-channel-joining
+    type: boolean?
+    default: true
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-squared-channel-joining'
+  - id: fit-spectral-pol
+    type: int?
+    default: 3
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-fit-spectral-pol'
+  - id: deconvolution-channels
+    type: int?
+    default: 3
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-deconvolution-channels'
+  - id: gridder
+    type: string?
+    default: wgridder
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-gridder'
+  - id: apply-primary-beam
+    type: boolean?
+    default: true
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-apply-primary-beam'
+  - id: use-differential-lofar-beam
+    type: boolean?
+    default: true
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-use-differential-lofar-beam'
+  - id: facet-regions
+    type: File?
+    inputBinding:
+      position: 1
+      shellQuote: false
+      prefix: '-facet-regions'
+
+  - id: facet-options
+    type:
+      type: record
+      name: facet_options
+      fields:
+        - name: facet-solutions
+          type: File?
+          inputBinding:
+            prefix: '-apply-facet-solutions'
+        - name: soltabs
+          type: string[]?
+          inputBinding:
+            itemSeparator: ','
+    default:
+      facet-solutions: null
+      soltabs: null
+
+outputs:
+  - id: Q_channel_images
+    type: File[]
+    doc: Per-channel Stokes Q images.
+    outputBinding:
+      glob: '$(inputs.name)-????-Q-image.fits'
+  - id: U_channel_images
+    type: File[]
+    doc: Per-channel Stokes U images.
+    outputBinding:
+      glob: '$(inputs.name)-????-U-image.fits'
+
+hints:
+  - class: DockerRequirement
+    dockerPull: vlbi-cwl
+
+requirements:
+  - class: ShellCommandRequirement
+  - class: InitialWorkDirRequirement
+    listing:
+      - entry: $(inputs.msin)
+  - class: ResourceRequirement
+    coresMin: $(inputs.cores)
+
+stdout: wsclean_qu.log
+stderr: wsclean_qu_err.log
diff --git a/workflows/polarization-imaging.cwl b/workflows/polarization-imaging.cwl
new file mode 100644
index 00000000..194fc3c9
--- /dev/null
+++ b/workflows/polarization-imaging.cwl
@@ -0,0 +1,155 @@
+class: Workflow
+cwlVersion: v1.2
+id: image_polarization 
+label: Polarization imaging
+doc: |
+  This workflow will image the provided MS in Q and U, and perform Rotation Measure Synthesis to provide linear polarization images.
+
+requirements:
+    - class: SubworkflowFeatureRequirement
+
+inputs:
+    - id: msin
+      type: Directory[]
+      doc: MeasurementSets that will be imaged.
+
+    - id: pixel_scale
+      type: string
+      doc: Pixel size (WSClean scale), e.g. "0.075asec".
+
+    - id: taper
+      type: string
+      doc: Angular resolution that will be passed to WSClean's taper argument. Its syntax follows that of WSClean.
+
+    - id: image_size
+      type: int[]
+      doc: Size (in pixels) of image, [x, y]. Its syntax follows that of WSClean.
+
+    - id: num_channels
+      type: int
+      doc: Number of channels to image in Q and U.
+
+    - id: stokes
+      type: string
+      default: "IQUV"
+
+    - id: rmtools_max_lam2
+      type: float?
+      default: 150
+      doc: Maximum lambda-squared value for rmsynth3d (-l).
+
+    - id: rmtools_dlam2
+      type: float?
+      default: 0.3
+      doc: Lambda-squared channel width for rmsynth3d (-d).
+
+    - id: rmtools_output_prefix
+      type: string?
+      doc: Prefix for RM-Tools output products. Defaults to Stokes Q basename.
+
+    - id: rmtools_extra_args
+      type: string?
+      doc: Extra arguments passed to rmsynth3d.
+
+steps:
+    - id: image_qu
+      label: image_qu
+      in:
+        - id: msin
+          source: msin
+        - id: scale
+          source: pixel_scale
+        - id: taper-gaussian
+          source: taper
+        - id: size
+          source: image_size
+        - id: channels-out
+          source: num_channels
+        - id: pol
+          source: stokes
+      out:
+        - id: Q_channel_images
+        - id: U_channel_images
+      run: ../steps/wsclean_pol.cwl
+
+    - id: make_cubes
+      label: Make QU cubes
+      in:
+        - id: Q_images
+          source: image_qu/Q_channel_images
+        - id: U_images
+          source: image_qu/U_channel_images
+        - id: image_size
+          source: image_size
+        - id: nchannels
+          source: num_channels
+      out:
+        - id: stokesQcube
+        - id: stokesUcube
+        - id: frequencies_list
+        - id: rms_list
+      run: ../steps/concat_pol.cwl
+
+    - id: run_rmtools
+      label: RM synthesis
+      in:
+        - id: stokes_q
+          source: make_cubes/stokesQcube
+        - id: stokes_u
+          source: make_cubes/stokesUcube
+        - id: freqs_hz
+          source: make_cubes/frequencies_list
+        - id: max_lam2
+          source: rmtools_max_lam2
+        - id: dlam2
+          source: rmtools_dlam2
+        - id: output_prefix
+          source: rmtools_output_prefix
+        - id: extra_args
+          source: rmtools_extra_args
+      out:
+        - id: fdf_im_dirty
+        - id: fdf_real_dirty
+        - id: fdf_tot_dirty
+        - id: fdf_maxpi
+        - id: fdf_peakrm
+        - id: rmsynth_stdout
+        - id: rmsynth_stderr
+      run: ../steps/run_rmtools.cwl
+
+outputs:
+  - id: stokesQcube
+    type: File
+    outputSource: make_cubes/stokesQcube
+
+  - id: stokesUcube
+    type: File
+    outputSource: make_cubes/stokesUcube
+
+  - id: fdf_im_dirty
+    type: File
+    outputSource: run_rmtools/fdf_im_dirty
+
+  - id: fdf_real_dirty
+    type: File
+    outputSource: run_rmtools/fdf_real_dirty
+
+  - id: fdf_tot_dirty
+    type: File
+    outputSource: run_rmtools/fdf_tot_dirty
+
+  - id: fdf_maxpi
+    type: File
+    outputSource: run_rmtools/fdf_maxpi
+
+  - id: fdf_peakrm
+    type: File
+    outputSource: run_rmtools/fdf_peakrm
+
+  - id: rmtools_stdout
+    type: File
+    outputSource: run_rmtools/rmsynth_stdout
+
+  - id: rmtools_stderr
+    type: File
+    outputSource: run_rmtools/rmsynth_stderr