diff --git a/examples/workflow/generic/multilevel_subworkflow_data/simple.ini b/examples/workflow/generic/multilevel_subworkflow_data/simple.ini
deleted file mode 120000
index 22f5e089deb..00000000000
--- a/examples/workflow/generic/multilevel_subworkflow_data/simple.ini
+++ /dev/null
@@ -1 +0,0 @@
-../simple_subworkflow_data/simple.ini
\ No newline at end of file
diff --git a/examples/workflow/generic/multilevel_subworkflow_data/simple.ini b/examples/workflow/generic/multilevel_subworkflow_data/simple.ini
new file mode 100644
index 00000000000..1cbe87c26c6
--- /dev/null
+++ b/examples/workflow/generic/multilevel_subworkflow_data/simple.ini
@@ -0,0 +1,9 @@
+[workflow]
+
+[executables]
+exe1 = /usr/bin/cp
+
+[exe1]
+
+[pegasus_profile]
+pycbc|submit-directory = ./
diff --git a/examples/workflow/generic/multilevel_subworkflow_data/simple.py b/examples/workflow/generic/multilevel_subworkflow_data/simple.py
deleted file mode 120000
index c8d9e43af52..00000000000
--- a/examples/workflow/generic/multilevel_subworkflow_data/simple.py
+++ /dev/null
@@ -1 +0,0 @@
-../simple_subworkflow_data/simple.py
\ No newline at end of file
diff --git a/examples/workflow/generic/multilevel_subworkflow_data/simple.py b/examples/workflow/generic/multilevel_subworkflow_data/simple.py
new file mode 100644
index 00000000000..f880258ed88
--- /dev/null
+++ b/examples/workflow/generic/multilevel_subworkflow_data/simple.py
@@ -0,0 +1,65 @@
+""" A minimal pycbc workflow example """
+
+import argparse
+import pycbc
+import pycbc.workflow as wf
+import os
+
+pycbc.init_logging(True)
+parser = argparse.ArgumentParser(description=__doc__[1:])
+parser.add_argument("--multilevel", action='store_true', default=False)
+wf.add_workflow_command_line_group(parser)
+wf.add_workflow_settings_cli(parser)
+args = parser.parse_args()
+
+input_file = wf.resolve_url_to_file("test_input.txt")
+input_file.add_pfn(os.path.abspath('./test_input.txt'), 'local')
+
+cont = wf.Workflow(args, 'cont')
+sub1 = wf.Workflow(args, 'sub1')
+sub1_1 = wf.Workflow(args, 'sub1_1')
+sub2 = wf.Workflow(args, 'sub2')
+
+exe1 = wf.Executable(cont.cp, 'exe1')
+
+SUBSUB = args.multilevel
+
+# Subworkflow 1: generate file that will be needed later
+
+# PATH1: generate that input in a sub-sub workflow
+if SUBSUB:
+    # Subworkflow 1: sub-subworkflow
+    node = exe1.create_node()
+    wf1_out_file = wf.File.from_path(os.path.abspath("test_output.txt.2"))
+    node.add_input_arg(input_file)
+    node.add_output_arg(wf1_out_file)
+    sub1_1 += node
+    sub1 += sub1_1
+    
+#PATH2: generate that input within a sub-workflow
+else:
+    node = exe1.create_node()
+    wf1_out_file = wf.File.from_path(os.path.abspath("test_output.txt.2"))
+    node.add_input_arg(input_file)
+    node.add_output_arg(wf1_out_file)
+    sub1 += node
+
+# Subworkflow 2
+# TEST GOAL: Job within subworkflow gets the input file produce in 
+# some external workflow
+node2 = exe1.create_node()
+node2.add_input_arg(wf1_out_file)
+node2.add_output_arg(wf.File.from_path(os.path.abspath("test_output.txt.3")))
+sub2 += node2
+
+# Regular job in top-level workflow
+# Test GOAL: job *directly* in workflow gets file produced by subworkflow in 
+# the same workflow
+node2 = exe1.create_node()
+node2.add_input_arg(wf1_out_file)
+node2.add_output_arg(wf.File.from_path(os.path.abspath("test_output.txt.4")))
+cont += node2
+
+cont += sub1
+cont += sub2
+cont.save()
diff --git a/examples/workflow/generic/multilevel_subworkflow_data/test_input.txt b/examples/workflow/generic/multilevel_subworkflow_data/test_input.txt
deleted file mode 120000
index 825b47e97fc..00000000000
--- a/examples/workflow/generic/multilevel_subworkflow_data/test_input.txt
+++ /dev/null
@@ -1 +0,0 @@
-../simple_subworkflow_data/test_input.txt
\ No newline at end of file
diff --git a/examples/workflow/generic/multilevel_subworkflow_data/test_input.txt b/examples/workflow/generic/multilevel_subworkflow_data/test_input.txt
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/examples/workflow/inference/small_test/single.ini b/examples/workflow/inference/small_test/single.ini
deleted file mode 120000
index b9e5ace0870..00000000000
--- a/examples/workflow/inference/small_test/single.ini
+++ /dev/null
@@ -1 +0,0 @@
-../../../inference/single/single_simple.ini
\ No newline at end of file
diff --git a/examples/workflow/inference/small_test/single.ini b/examples/workflow/inference/small_test/single.ini
new file mode 100644
index 00000000000..16d71b95651
--- /dev/null
+++ b/examples/workflow/inference/small_test/single.ini
@@ -0,0 +1,66 @@
+[model]
+name = single_template
+
+#; This model precalculates the SNR time series at a fixed rate.
+#; If you need a higher time resolution, this may be increased
+sample_rate = 32768
+low-frequency-cutoff = 30.0
+
+[data]
+instruments = H1 L1 V1
+analysis-start-time = 1187008482
+analysis-end-time = 1187008892
+psd-estimation = median
+psd-segment-length = 16
+psd-segment-stride = 8
+psd-inverse-length = 16
+pad-data = 8
+channel-name = H1:LOSC-STRAIN L1:LOSC-STRAIN V1:LOSC-STRAIN
+frame-files = H1:H-H1_LOSC_CLN_4_V1-1187007040-2048.gwf L1:L-L1_LOSC_CLN_4_V1-1187007040-2048.gwf V1:V-V1_LOSC_CLN_4_V1-1187007040-2048.gwf
+strain-high-pass = 15
+sample-rate = 2048
+
+[sampler]
+name = dynesty
+sample = rwalk
+bound = multi
+dlogz = 0.1
+nlive = 200
+checkpoint_time_interval = 100
+maxcall = 10000
+
+[variable_params]
+; waveform parameters that will vary in MCMC
+tc =
+distance =
+inclination =
+
+[static_params]
+; waveform parameters that will not change in MCMC
+approximant = TaylorF2
+f_lower = 30
+mass1 = 1.3757
+mass2 = 1.3757
+
+#; we'll choose not to sample over these, but you could
+polarization = 0
+ra = 3.44615914
+dec = -0.40808407
+
+#; You could also set additional parameters if your waveform model supports / requires it.
+; spin1z = 0
+
+[prior-tc]
+; coalescence time prior
+name = uniform
+min-tc = 1187008882.4
+max-tc = 1187008882.5
+
+[prior-distance]
+#; following gives a uniform in volume
+name = uniform_radius
+min-distance = 10
+max-distance = 60
+
+[prior-inclination]
+name = sin_angle
diff --git a/pycbc/conversions.py b/pycbc/conversions.py
index 04f5dc1779d..311eaafb05b 100644
--- a/pycbc/conversions.py
+++ b/pycbc/conversions.py
@@ -1198,6 +1198,47 @@ def get_lm_f0tau_allmodes(mass, spin, modes):
             tau[key.format(l, abs(m), n)] = tmp_tau
     return f0, tau
 
+def get_qlm_f0tau_allmodes(mass, spin, qmodes):
+    """Returns a dictionary of all of the frequencies and damping times for the
+    requested modes.
+
+    Parameters
+    ----------
+    mass : float or array
+        Mass of the black hole (in solar masses).
+    spin : float or array
+        Dimensionless spin of the final black hole.
+    qmodes : list of str
+        The modes to get. Each string in the list should be formatted
+        'qlmN', where l (m) is the l (m) index of the
+        harmonic and N is the number of overtones to generate (note, N is not
+        the index of the overtone).
+
+    Returns
+    -------
+    f0 : dict
+        Dictionary mapping the modes to the frequencies. The dictionary keys
+        are 'qlmn' string, where l (m) is the l (m) index of the harmonic and
+        n is the index of the overtone. For example, '220' is the l = m = 2
+        mode and the 0th overtone.
+    tau : dict
+        Dictionary mapping the modes to the damping times. The keys are the
+        same as ``f0``.
+    """
+    f0, tau = {}, {}
+    for qlmn, lmn1, lmn2 in qmodes:
+        key = '{}{}{}'
+        l1, m1, nmodes1 = int(lmn1[0]), int(lmn1[1]), int(lmn1[2])
+        l2, m2, nmodes2 = int(lmn2[0]), int(lmn2[1]), int(lmn2[2])
+        ql, qm, qnmodes = int(qlmn[0]), int(qlmn[1]), int(qlmn[2])
+        for qn in range(qnmodes):
+            for n1 in range(nmodes1):
+                for n2 in range(nmodes2):
+                    tmp_f01, tmp_tau1 = get_lm_f0tau(mass, spin, l1, m1, n1)
+                    tmp_f02, tmp_tau2 = get_lm_f0tau(mass, spin, l2, m2, n2)
+                    f0[key.format(ql, abs(qm), qn)] = tmp_f01 + tmp_f02
+                    tau[key.format(ql, abs(qm), qn)] = (tmp_tau1 * tmp_tau2) / (tmp_tau1 + tmp_tau2)
+    return f0, tau
 
 def freq_from_final_mass_spin(final_mass, final_spin, l=2, m=2, n=0):
     """Returns QNM frequency for the given mass and spin and mode.
diff --git a/pycbc/waveform/ringdown.py b/pycbc/waveform/ringdown.py
index 6f4feca6dab..2a1512c33fb 100644
--- a/pycbc/waveform/ringdown.py
+++ b/pycbc/waveform/ringdown.py
@@ -33,11 +33,11 @@
 from pycbc.types import (TimeSeries, FrequencySeries, float64, complex128,
                          zeros)
 from pycbc.waveform.waveform import get_obj_attrs
-from pycbc.conversions import get_lm_f0tau_allmodes
+from pycbc.conversions import get_lm_f0tau_allmodes, get_qlm_f0tau_allmodes
 
 qnm_required_args = ['f_0', 'tau', 'amp', 'phi']
-mass_spin_required_args = ['final_mass','final_spin', 'lmns', 'inclination']
-freqtau_required_args = ['lmns']
+mass_spin_required_args = ['final_mass','final_spin', 'inclination']
+freqtau_required_args = []
 td_args = {'delta_t': None, 't_final': None, 'taper': False}
 fd_args = {'t_0': 0, 'delta_f': None, 'f_lower': 0, 'f_final': None}
 
@@ -119,6 +119,50 @@ def format_lmns(lmns):
 
     return out
 
+def format_qlmns(qlmns):
+    """Checks if the format of the parameter qlmns is correct, returning the
+    appropriate format if not, and raise an error if nmodes=0.
+
+    The required format for the ringdown approximants is a list of qlmn modes
+    as a single whitespace-separated string, with n the number
+    of overtones desired. Alternatively, a list object may be used,
+    containing individual strings as elements for the qlmn modes.
+    For instance, qlmns = '221 331' are the modes 220, 221, 222, and 330.
+    Giving qlmns = ['221', '331'] is equivalent.
+    The ConfigParser of a workflow might convert that to a single string
+    (case 1 below) or a list with a single string (case 2), and this function
+    will return the appropriate list of strings. If a different format is
+    given, raise an error.
+    """
+
+    # Catch case of qlmns given as float (as int injection values are cast
+    # to float by pycbc_create_injections), cast to int, then string
+    if isinstance(qlmns, float):
+        qlmns = str(int(qlmns))
+    # Case 1: the qlmns are given as a string, e.g. '221 331'
+    if isinstance(qlmns, str):
+        qlmns = qlmns.split(' ')
+    # Case 2: the qlmns are given as strings in a list, e.g. ['221', '331']
+    elif isinstance(qlmns, list):
+        pass
+    else:
+        raise ValueError('Format of parameter qlmns not recognized. See '
+                         'approximant documentation for more info.')
+
+    out = []
+    # Cycle over the qlmns to ensure that we get back a list of strings that
+    # are three digits long, and that nmodes!=0
+    for qlmn in qlmns:
+        # Quadratic modes are given as qlmn=lmn1xlmn2
+        qlmn, lmn12 = qlmn.split('=')
+        lmn1, lmn2 = lmn12.split('x')
+        if not (abs(int(lmn1[0])-int(lmn2[0]))<= int(qlmn[0]) <= int(lmn1[0])+int(lmn2[0])):
+            raise ValueError('Quadratic l={} must be in the range |l1-l2|<= l <= l1+l2'.format(qlmn[0]))
+        if not (abs(int(lmn1[1]))+abs(int(lmn2[1])) == abs(int(qlmn[1]))):
+            raise ValueError('Quadratic m={} must satisfy |m1|+|m2| = |m|'.format(qlmn[1]))
+        out.append((qlmn, lmn1, lmn2))
+    return out
+
 def parse_mode(lmn):
     """Extracts overtones from an lmn.
     """
@@ -179,6 +223,70 @@ def lm_amps_phases(**kwargs):
             dbetas[mode] = kwargs.pop('dbeta'+mode, ref_dbeta)
     return amps, phis, dbetas, dphis
 
+def qlm_amps_phases(**kwargs):
+    r"""Takes input_params and return dictionaries with amplitudes and phases
+    of each overtone of a specific qlm mode, checking that all of them are
+    given. Will also look for dbetas and dphis. If ``(dphi|dbeta)`` (i.e.,
+    without a mode suffix) are provided, they will be used for all modes that
+    don't explicitly set a ``(dphi|dbeta){lmn}``.
+    """
+    amps = {}
+    phis = {}
+    dbetas = {}
+    dphis = {}
+    qref_amp = kwargs.pop('qref_amp', None)
+    if qref_amp is None:
+        # default to the 220 mode
+        qref_amp = 1
+    qref_phi = kwargs.pop('qref_phi', None)
+    if qref_phi is None:
+        # default to the 220 mode
+        qref_phi = 0
+    # reference mode
+    ref_amp = kwargs.pop('ref_amp', None)
+    if ref_amp is None:
+        # default to the 220 mode
+        ref_amp = 'amp220'
+    # check for reference dphi and dbeta
+    ref_dbeta = kwargs.pop('dbeta', 0.)
+    ref_dphi = kwargs.pop('dphi', 0.)
+    if isinstance(ref_amp, str) and ref_amp.startswith('amp'):
+        # assume a mode was provided; check if the mode exists
+        ref_mode = ref_amp.replace('amp', '')
+        try:
+            ref_amp = kwargs.pop(ref_amp)
+            amps[ref_mode] = ref_amp
+        except KeyError:
+            raise ValueError("Must provide an amplitude for the reference "
+                             "mode {}".format(ref_amp))
+    else:
+        ref_mode = None
+    # Get amplitudes and phases of the modes
+    for qlmn, lmn1, lmn2 in kwargs['qlmns']:
+        overtones_qlmn = parse_mode(qlmn)
+        overtones_lmn1, overtones_lmn2 = parse_mode(lmn1), parse_mode(lmn2)
+        for modeq, mode1, mode2 in zip(overtones_qlmn, overtones_lmn1, overtones_lmn2):
+            try:
+                if mode1 != ref_mode and mode2 != ref_mode:
+                    amps[modeq] = qref_amp * kwargs['amp' + mode1] * kwargs['amp' + mode2]
+                elif mode1 == ref_mode and mode2 != ref_mode:
+                    amps[modeq] = qref_amp * ref_amp * kwargs['amp' + mode2]
+                elif mode1 != ref_mode and mode2 == ref_mode:
+                    amps[modeq] = qref_amp * kwargs['amp' + mode1] * ref_amp
+                else:
+                    amps[modeq] = qref_amp * ref_amp * ref_amp
+            except KeyError:
+                raise ValueError('Must provide reference amplitude for quadratic modes, and \
+                                    amp{} and amp{} are required'.format(mode1, mode2))
+            try:
+                phis[modeq] = kwargs['phi' + mode1] + kwargs['phi' + mode2] + qref_phi
+            except KeyError:
+                raise ValueError('Must provide reference phi for quadratic modes, and \
+                                    phi{} and phi{} are required'.format(mode1, mode2))
+            dphis[modeq] = kwargs.pop('dphi'+qlmn, ref_dphi)
+            dbetas[modeq] = kwargs.pop('dbeta'+qlmn, ref_dbeta)                
+    return amps, phis, dbetas, dphis
+
 
 def lm_freqs_taus(**kwargs):
     """Take input_params and return dictionaries with frequencies and damping
@@ -200,6 +308,27 @@ def lm_freqs_taus(**kwargs):
                 raise ValueError('tau_{} is required'.format(mode))
     return freqs, taus
 
+def qlm_freqs_taus(**kwargs):
+    """Take input_params and return dictionaries with frequencies and damping
+    times of each overtone of a specific qlm mode, checking that all of them
+    are given.
+    """
+    freqs, taus = {}, {}
+    for qlmn, lmn1, lmn2 in kwargs['qlmns']:
+        overtones_qlmn = parse_mode(qlmn)
+        overtones_lmn1, overtones_lmn2 = parse_mode(lmn1), parse_mode(lmn2)
+        for modeq, mode1, mode2 in zip(overtones_qlmn, overtones_lmn1, overtones_lmn2):
+            try:
+                freqs[modeq] = kwargs['f_' + mode1] + kwargs['f_' + mode2]
+            except KeyError:
+                raise ValueError('f_{} and f_{} are required'.format(mode1, mode2))
+            try:
+                taus[modeq] = (kwargs['tau_' + mode1] * kwargs['tau_' + mode2]) / \
+                              (kwargs['tau_' + mode1] + kwargs['tau_' + mode2])
+            except KeyError:
+                raise ValueError('tau_{} and tau_{} are required'.format(mode1, mode2))
+    return freqs, taus
+
 
 def lm_arbitrary_harmonics(**kwargs):
     """Take input_params and return dictionaries with arbitrary harmonics
@@ -215,6 +344,20 @@ def lm_arbitrary_harmonics(**kwargs):
             polnms[mode] = kwargs.pop('polnm{}'.format(mode), None)
     return pols, polnms
 
+def qlm_arbitrary_harmonics(**kwargs):
+    """Take input_params and return dictionaries with arbitrary harmonics
+    for each mode.
+    """
+    pols = {}
+    polnms = {}
+    for qlmn, lmn1, lmn2 in kwargs['qlmns']:
+        overtones_qlmn = parse_mode(qlmn)
+        overtones_lmn1, overtones_lmn2 = parse_mode(lmn1), parse_mode(lmn2)
+        for modeq, mode1, mode2 in zip(overtones_qlmn, overtones_lmn1, overtones_lmn2):
+            pols[modeq] = kwargs.pop('pol{}'.format(mode1), None)
+            polnms[modeq] = kwargs.pop('polnm{}'.format(mode2), None)
+    return pols, polnms
+
 
 # Functions to obtain t_final, f_final and output vector ######################
 
@@ -729,9 +872,18 @@ def multimode_base(input_params, domain, freq_tau_approximant=False):
         The cross phase of a ringdown with the lm modes specified and
         n overtones in the chosen domain (time or frequency).
     """
-    input_params['lmns'] = format_lmns(input_params['lmns'])
-    amps, phis, dbetas, dphis = lm_amps_phases(**input_params)
-    pols, polnms = lm_arbitrary_harmonics(**input_params)
+    try:    
+        input_params['lmns'] = format_lmns(input_params['lmns'])
+        amps, phis, dbetas, dphis = lm_amps_phases(**input_params)
+        pols, polnms = lm_arbitrary_harmonics(**input_params)
+    except KeyError:
+        input_params['lmns'] = None
+    try:    
+        input_params['qlmns'] = format_qlmns(input_params['qlmns'])
+        q_amps, q_phis, q_dbetas, q_dphis = qlm_amps_phases(**input_params)
+        q_pols, q_polnms = qlm_arbitrary_harmonics(**input_params)
+    except KeyError:
+        input_params['qlmns'] = None    
     # get harmonics argument
     try:
         harmonics = input_params['harmonics']
@@ -749,63 +901,130 @@ def multimode_base(input_params, domain, freq_tau_approximant=False):
         input_params['azimuthal'] = 0.
     # figure out the frequencies and damping times
     if freq_tau_approximant:
-        freqs, taus = lm_freqs_taus(**input_params)
+        if input_params['lmns'] is not None:
+            freqs, taus = lm_freqs_taus(**input_params)
+        if input_params['qlmns'] is not None:
+            qfreqs, qtaus = qlm_freqs_taus(**input_params)
         norm = 1.
     else:
-        freqs, taus = get_lm_f0tau_allmodes(input_params['final_mass'],
-                        input_params['final_spin'], input_params['lmns'])
-        norm = Kerr_factor(input_params['final_mass'],
-            input_params['distance']) if 'distance' in input_params.keys() \
-            else 1.
-        for mode, freq in freqs.items():
-            if 'delta_f{}'.format(mode) in input_params:
-                freqs[mode] += input_params['delta_f{}'.format(mode)]*freq
-        for mode, tau in taus.items():
-            if 'delta_tau{}'.format(mode) in input_params:
-                taus[mode] += input_params['delta_tau{}'.format(mode)]*tau
+        if input_params['lmns'] is not None:
+            freqs, taus = get_lm_f0tau_allmodes(input_params['final_mass'],
+                            input_params['final_spin'], input_params['lmns'])            
+            norm = Kerr_factor(input_params['final_mass'],
+                input_params['distance']) if 'distance' in input_params.keys() \
+                else 1.
+            for mode, freq in freqs.items():
+                if 'delta_f{}'.format(mode) in input_params:
+                    freqs[mode] += input_params['delta_f{}'.format(mode)]*freq
+            for mode, tau in taus.items():
+                if 'delta_tau{}'.format(mode) in input_params:
+                    taus[mode] += input_params['delta_tau{}'.format(mode)]*tau
+        if input_params['qlmns'] is not None:
+            qfreqs, qtaus = get_qlm_f0tau_allmodes(input_params['final_mass'],
+                            input_params['final_spin'], input_params['qlmns'])
+            norm = Kerr_factor(input_params['final_mass'],
+                input_params['distance']) if 'distance' in input_params.keys() \
+                else 1.
+            for mode, qfreq in qfreqs.items():
+                if 'delta_f{}'.format(mode) in input_params:
+                    qfreqs[mode] += input_params['delta_f{}'.format(mode)]*qfreq
+            for mode, qtau in qtaus.items():
+                if 'delta_tau{}'.format(mode) in input_params:
+                    qtaus[mode] += input_params['delta_tau{}'.format(mode)]*qtau
     # setup the output
     if domain == 'td':
-        outplus, outcross = td_output_vector(freqs, taus,
+        if input_params['lmns'] is not None:
+            if input_params['qlmns'] is not None:
+                outplus, outcross = td_output_vector(qfreqs, qtaus,
+                                input_params['taper'], input_params['delta_t'],
+                                input_params['t_final'])
+            else:
+                outplus, outcross = td_output_vector(freqs, taus,
+                                input_params['taper'], input_params['delta_t'],
+                                input_params['t_final'])
+            sample_times = outplus.sample_times.numpy()
+        elif input_params['qlmns'] is not None:
+            outplus, outcross = td_output_vector(qfreqs, qtaus,
                             input_params['taper'], input_params['delta_t'],
                             input_params['t_final'])
-        sample_times = outplus.sample_times.numpy()
+            sample_times = outplus.sample_times.numpy()
     elif domain == 'fd':
-        outplus, outcross = fd_output_vector(freqs, taus,
+        if input_params['lmns'] is not None:
+            if input_params['qlmns'] is not None:
+                outplus, outcross = fd_output_vector(qfreqs, qtaus,
+                                input_params['delta_f'], input_params['f_final'])
+            else:
+                outplus, outcross = fd_output_vector(freqs, taus,
+                                input_params['delta_f'], input_params['f_final'])
+            kmin = int(input_params['f_lower'] / input_params['delta_f'])
+            sample_freqs = outplus.sample_frequencies.numpy()[kmin:]
+        elif input_params['qlmns'] is not None:
+            outplus, outcross = fd_output_vector(qfreqs, qtaus,
                             input_params['delta_f'], input_params['f_final'])
-        kmin = int(input_params['f_lower'] / input_params['delta_f'])
-        sample_freqs = outplus.sample_frequencies.numpy()[kmin:]
+            kmin = int(input_params['f_lower'] / input_params['delta_f'])
+            sample_freqs = outplus.sample_frequencies.numpy()[kmin:]
     else:
         raise ValueError('unrecognised domain argument {}; '
                          'must be either fd or td'.format(domain))
     # cyclce over the modes, generating the waveforms
-    for lmn in freqs:
-        if amps[lmn] == 0.:
-            # skip
-            continue
-        if domain == 'td':
-            hplus, hcross = td_damped_sinusoid(
-                freqs[lmn], taus[lmn], amps[lmn], phis[lmn], sample_times,
-                l=int(lmn[0]), m=int(lmn[1]), n=int(lmn[2]),
-                inclination=input_params['inclination'],
-                azimuthal=input_params['azimuthal'],
-                dphi=dphis[lmn], dbeta=dbetas[lmn],
-                harmonics=harmonics, final_spin=final_spin,
-                pol=pols[lmn], polnm=polnms[lmn])
-            outplus += hplus
-            outcross += hcross
-        elif domain == 'fd':
-            hplus, hcross = fd_damped_sinusoid(
-                freqs[lmn], taus[lmn], amps[lmn], phis[lmn], sample_freqs,
-                l=int(lmn[0]), m=int(lmn[1]), n=int(lmn[2]),
-                inclination=input_params['inclination'],
-                azimuthal=input_params['azimuthal'],
-                harmonics=harmonics, final_spin=final_spin,
-                pol=pols[lmn], polnm=polnms[lmn])
-            outplus[kmin:] += hplus
-            outcross[kmin:] += hcross
+    if input_params['lmns'] is not None:
+        for lmn in freqs:
+            if amps[lmn] == 0.:
+                # skip
+                continue
+            if domain == 'td':
+                print('{}: Frequency: {}, Tau: {}'.format(lmn, freqs[lmn], taus[lmn]))
+                hplus, hcross = td_damped_sinusoid(
+                    freqs[lmn], taus[lmn], amps[lmn], phis[lmn], sample_times,
+                    l=int(lmn[0]), m=int(lmn[1]), n=int(lmn[2]),
+                    inclination=input_params['inclination'],
+                    azimuthal=input_params['azimuthal'],
+                    dphi=dphis[lmn], dbeta=dbetas[lmn],
+                    harmonics=harmonics, final_spin=final_spin,
+                    pol=pols[lmn], polnm=polnms[lmn])
+                outplus += hplus
+                outcross += hcross
+            elif domain == 'fd':
+                print('{}: Frequency: {}, Tau: {}'.format(lmn, freqs[lmn], taus[lmn]))
+                hplus, hcross = fd_damped_sinusoid(
+                    freqs[lmn], taus[lmn], amps[lmn], phis[lmn], sample_freqs,
+                    l=int(lmn[0]), m=int(lmn[1]), n=int(lmn[2]),
+                    inclination=input_params['inclination'],
+                    azimuthal=input_params['azimuthal'],
+                    harmonics=harmonics, final_spin=final_spin,
+                    pol=pols[lmn], polnm=polnms[lmn])
+                outplus[kmin:] += hplus
+                outcross[kmin:] += hcross
+    if input_params['qlmns'] is not None:
+        for qlmn in qfreqs:
+            if q_amps[qlmn] == 0.:
+                # skip
+                continue
+            if domain == 'td':
+                print('{}: Frequency: {}, Tau: {}'.format(qlmn, qfreqs[qlmn], qtaus[qlmn]))
+                hplus, hcross = td_damped_sinusoid(
+                    qfreqs[qlmn], qtaus[qlmn], q_amps[qlmn], q_phis[qlmn],
+                    sample_times, l=int(qlmn[0]), m=int(qlmn[1]), n=int(qlmn[2]),
+                    inclination=input_params['inclination'],
+                    azimuthal=input_params['azimuthal'],
+                    dphi=q_dphis[qlmn], dbeta=q_dbetas[qlmn],
+                    harmonics=harmonics, final_spin=final_spin,
+                    pol=q_pols[qlmn], polnm=q_polnms[qlmn])
+                outplus += hplus
+                outcross += hcross
+            elif domain == 'fd':
+                print('{}: Frequency: {}, Tau: {}'.format(qlmn, qfreqs[qlmn], qtaus[qlmn]))
+                hplus, hcross = fd_damped_sinusoid(
+                    qfreqs[qlmn], qtaus[qlmn], q_amps[qlmn], q_phis[qlmn],
+                    sample_freqs, l=int(qlmn[0]), m=int(qlmn[1]), n=int(qlmn[2]),
+                    inclination=input_params['inclination'],
+                    azimuthal=input_params['azimuthal'],
+                    harmonics=harmonics, final_spin=final_spin,
+                    pol=q_pols[qlmn], polnm=q_polnms[qlmn])
+                outplus[kmin:] += hplus
+                outcross[kmin:] += hcross
     return norm * outplus, norm * outcross
 
-
 ######################################################
 #### Approximants
 ######################################################
@@ -1231,4 +1450,4 @@ def get_fd_from_freqtau(template=None, **kwargs):
 
 ringdown_td_approximants = {
     'TdQNMfromFinalMassSpin': get_td_from_final_mass_spin,
-    'TdQNMfromFreqTau': get_td_from_freqtau}
+    'TdQNMfromFreqTau': get_td_from_freqtau}
\ No newline at end of file