diff --git a/sotodlib/mapmaking/demod_mapmaker.py b/sotodlib/mapmaking/demod_mapmaker.py
index 8b13d9152..0eebc7d61 100644
--- a/sotodlib/mapmaking/demod_mapmaker.py
+++ b/sotodlib/mapmaking/demod_mapmaker.py
@@ -527,7 +527,7 @@ def make_demod_map(context, obslist, noise_model, info,
                                                 overwrite=False)
         errors.append(error) ; outputs.append((output_init, output_proc)) ;
         if error not in [None,'load_success']:
-            L.info('tod %s:%s:%s failed in the prepoc database'%(obs_id,detset,band))
+            L.info('tod %s:%s:%s failed in the preproc database'%(obs_id,detset,band))
             continue
         obs.wrap("weather", np.full(1, "toco"))
         obs.wrap("site",    np.full(1, site))
diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py
index 488b18a57..2d407b507 100644
--- a/sotodlib/preprocess/processes.py
+++ b/sotodlib/preprocess/processes.py
@@ -804,7 +804,7 @@ def process(self, aman, proc_aman):
                                     aman.samps.offset + aman.samps.count - trim))
 
 
-class EstimateAzSS(_Preprocess):
+class AzSS(_Preprocess):
     """Estimates Azimuth Synchronous Signal (AzSS) by binning signal by azimuth of boresight.
     All process confgis go to `get_azss`. If `method` is 'interpolate', no fitting applied 
     and binned signal is directly used as AzSS model. If `method` is 'fit', Legendre polynominal
@@ -812,29 +812,68 @@ class EstimateAzSS(_Preprocess):
 
     Example configuration block::
 
-      - name: "estimate_azss"
+      - name: "azss"
+        azss_stats_name: 'azss_statsQ'
+        proc_aman_turnaround_info: 'turnaround_flags'
         calc:
           signal: 'demodQ'
-          azss_stats_name: 'azss_statsQ'
-          range: [-1.57079, 7.85398]
+          frange: [-1.57079, 7.85398]
           bins: 1080
-          merge_stats: False
-          merge_model: False
+          left_right: True
         save: True
+        process:
+            subtract: True
 
     .. autofunction:: sotodlib.tod_ops.azss.get_azss
     """
-    name = "estimate_azss"
+    name = "azss"
+    def __init__(self, step_cfgs):
+        self.azss_stats_name = step_cfgs.get('azss_stats_name', 'azss_stats')
+        self.proc_aman_turnaround_info = step_cfgs.get('proc_aman_turnaround_info', None)
+
+        super().__init__(step_cfgs)
 
     def calc_and_save(self, aman, proc_aman):
-        calc_aman, _ = tod_ops.azss.get_azss(aman, **self.calc_cfgs)
-        self.save(proc_aman, calc_aman)
+        # If process is run then just wrap info from process step
+        if self.process_cfgs:
+            self.save(proc_aman, aman[self.azss_stats_name])
+        else:
+            if self.proc_aman_turnaround_info:
+                _f = attrgetter(self.proc_aman_turnaround_info)
+                turnaround_info = _f(proc_aman)
+            else:
+                turnaround_info = None
+            _azss = tod_ops.azss.get_azss(aman, azss_stats_name=self.azss_stats_name,
+                                  turnaround_info=turnaround_info,
+                                  merge_stats=False, merge_model=False,
+                                  subtract_in_place=self.process_cfgs["subtract"], 
+                                  **self.calc_cfgs)
+            azss_stats = _azss[0]
+            self.save(proc_aman, azss_stats)
     
     def save(self, proc_aman, azss_stats):
         if self.save_cfgs is None:
             return
         if self.save_cfgs:
-            proc_aman.wrap(self.calc_cfgs["azss_stats_name"], azss_stats)
+            proc_aman.wrap(self.azss_stats_name, azss_stats)
+
+    def process(self, aman, proc_aman):
+        if self.proc_aman_turnaround_info:
+            _f = attrgetter(self.proc_aman_turnaround_info)
+            turnaround_info = _f(proc_aman)
+        else:
+            turnaround_info = None
+        if self.azss_stats_name in proc_aman:
+            if self.process_cfgs["subtract"]:
+                tod_ops.azss.subtract_azss(aman, proc_aman[self.azss_stats_name],
+                                           signal = self.calc_cfgs.get('signal'),
+                                           in_place=True)
+        else:
+            tod_ops.azss.get_azss(aman, azss_stats_name=self.azss_stats_name,
+                                  turnaround_info=turnaround_info,
+                                  merge_stats=True, merge_model=False,
+                                  subtract_in_place=self.process_cfgs["subtract"], 
+                                  **self.calc_cfgs)
 
 class GlitchFill(_Preprocess):
     """Fill glitches. All process configs go to `fill_glitches`.
@@ -1546,6 +1585,7 @@ def __init__(self, step_cfgs):
         super().__init__(step_cfgs)
 
     def calc_and_save(self, aman, proc_aman):
+        coeff_aman = get_qu_common_mode_coeffs(aman, Q_signal, U_signal, merge)
         self.save(proc_aman, aman)
 
     def save(self, proc_aman, aman):
@@ -1592,7 +1632,7 @@ def save(self, proc_aman, fp_aman):
             return
         if self.save_cfgs:
             proc_aman.wrap("fp_flags", fp_aman)
-    
+
     def select(self, meta, proc_aman=None):
         if self.select_cfgs is None:
             return meta
@@ -1621,6 +1661,60 @@ def process(self, aman, proc_aman):
         if self.process_cfgs:
             pointing_model.apply_pointing_model(aman)
 
+class BadSubscanFlags(_Preprocess):
+    """Identifies and flags bad subscans.
+
+      Example config block::
+
+        - name : "noisy_subscan_flags"
+            stats_name: tod_stats # optional
+            calc: 
+                nstd_lim: 5.0 
+                merge: False
+            save: True
+            select: True
+    
+    .. autofunction:: sotodlib.tod_ops.flags.get_badsubscan_flags
+    """
+    name = "noisy_subscan_flags"
+
+    def __init__(self, step_cfgs):
+        self.stats_name = step_cfgs.get('stats_name', 'tod_stats')
+        super().__init__(step_cfgs)
+    
+    def calc_and_save(self, aman, proc_aman):
+        from so3g.proj import RangesMatrix, Ranges
+        msk_ss = RangesMatrix.zeros((aman.dets.count, aman.samps.count))
+        msk_det = np.ones(aman.dets.count, dtype=bool)
+        signal = self.calc_cfgs["subscan_stats"] if isinstance(self.calc_cfgs["subscan_stats"], list) else [self.calc_cfgs["subscan_stats"]]
+        for sig in signal:
+            self.calc_cfgs["subscan_stats"] = proc_aman[self.stats_name+"_"+sig[-1]]
+            _msk_ss, _msk_det = tod_ops.flags.get_noisy_subscan_flags(
+                aman, **self.calc_cfgs)
+            msk_ss += _msk_ss
+            msk_det &= _msk_det
+        ss_aman = core.AxisManager(aman.dets, aman.samps)
+        ss_aman.wrap("valid_subscans", msk_ss, [(0, 'dets'), (1, 'samps')])
+        det_aman = core.AxisManager(aman.dets)
+        det_aman.wrap("valid_dets", msk_det)
+        self.save(proc_aman, ss_aman, "noisy_subscan_flags")
+        self.save(proc_aman, det_aman, "noisy_dets_flags")
+
+    def save(self, proc_aman, calc_aman, name): 
+        if self.save_cfgs is None:
+            return
+        if self.save_cfgs:
+            proc_aman.wrap(name, calc_aman)
+
+    def select(self, meta, proc_aman=None):
+        if self.select_cfgs is None:
+            return meta
+        if proc_aman is None:
+            proc_aman = meta.preprocess
+        keep = proc_aman.noisy_dets_flags.valid_dets
+        meta.restrict('dets', proc_aman.dets.vals[keep])
+        return meta
+
 _Preprocess.register(SplitFlags)
 _Preprocess.register(SubtractT2P)
 _Preprocess.register(EstimateT2P)
@@ -1641,7 +1735,7 @@ def process(self, aman, proc_aman):
 _Preprocess.register(SubtractHWPSS)
 _Preprocess.register(Apodize)
 _Preprocess.register(Demodulate)
-_Preprocess.register(EstimateAzSS)
+_Preprocess.register(AzSS)
 _Preprocess.register(GlitchFill)
 _Preprocess.register(FlagTurnarounds)
 _Preprocess.register(SubPolyf)
@@ -1655,4 +1749,5 @@ def process(self, aman, proc_aman):
 _Preprocess.register(RotateQU) 
 _Preprocess.register(SubtractQUCommonMode) 
 _Preprocess.register(FocalplaneNanFlags) 
-_Preprocess.register(PointingModel) 
+_Preprocess.register(PointingModel)  
+_Preprocess.register(BadSubscanFlags)
diff --git a/sotodlib/tod_ops/azss.py b/sotodlib/tod_ops/azss.py
index e95f3ad43..4993939cc 100644
--- a/sotodlib/tod_ops/azss.py
+++ b/sotodlib/tod_ops/azss.py
@@ -1,15 +1,17 @@
 """Module for estimating Azimuth Synchronous Signal (azss)"""
 import numpy as np
+from operator import attrgetter
 from numpy.polynomial import legendre as L
 from scipy.optimize import curve_fit
 from scipy.interpolate import interp1d
 from sotodlib import core, tod_ops
 from sotodlib.tod_ops import bin_signal, apodize, filters
+from so3g.proj import Ranges, RangesMatrix
 import logging
 
 logger = logging.getLogger(__name__)
 
-def bin_by_az(aman, signal=None, az=None, range=None, bins=100, flags=None, 
+def bin_by_az(aman, signal=None, az=None, frange=None, bins=100, flags=None, 
               apodize_edges=True, apodize_edges_samps=1600, 
               apodize_flags=True, apodize_flags_samps=200):
     """
@@ -23,13 +25,13 @@ def bin_by_az(aman, signal=None, az=None, range=None, bins=100, flags=None,
         numpy array of signal to be binned. If None, the signal is taken from aman.signal.
     az: array-like, optional
         A 1D numpy array representing the azimuth angles. If not provided, the azimuth angles are taken from aman.boresight.az attribute.
-    range: array-like, optional
+    frange: array-like, optional
         A list specifying the range of azimuth angles to consider for binning. Defaults to None.
         If None, [min(az), max(az)] will be used for binning.
     bins: int or sequence of scalars
         If bins is an int, it defines the number of equal-width bins in the given range (100, by default).
         If bins is a sequence, it defines the bin edges, including the rightmost edge, allowing for non-uniform bin widths.
-        If ``bins`` is a sequence, ``bins`` overwrite ``range``.
+        If ``bins`` is a sequence, ``bins`` overwrite ``frange``.
     flags: RangesMatrix, optional
         Flag indicating whether to exclude flagged samples when binning the signal.
         Default is no mask applied.
@@ -80,7 +82,7 @@ def bin_by_az(aman, signal=None, az=None, range=None, bins=100, flags=None,
         else:
             weight_for_signal = None
     binning_dict = bin_signal(aman, bin_by=az, signal=signal,
-                               range=range, bins=bins, flags=flags, weight_for_signal=weight_for_signal)
+                              range=frange, bins=bins, flags=flags, weight_for_signal=weight_for_signal)
     return binning_dict
 
 def fit_azss(az, azss_stats, max_mode, fit_range=None):
@@ -141,13 +143,42 @@ def fit_azss(az, azss_stats, max_mode, fit_range=None):
     
     return azss_stats, L.legval(x_legendre, coeffs.T)
     
-    
-def get_azss(aman, signal='signal', az=None, range=None, bins=100, flags=None,
-            apodize_edges=True, apodize_edges_samps=40000, apodize_flags=True, apodize_flags_samps=200,
+def _prepare_azss_stats(aman, signal, az, frange=None, bins=100, flags=None, 
+                        apodize_edges=True, apodize_edges_samps=40000, apodize_flags=True,
+                        apodize_flags_samps=200, method='interpolate', max_mode=None):
+    """
+    Helper function to collect initial info for azss_stats AxisManager.
+    """
+    # do binning
+    binning_dict = bin_by_az(aman, signal=signal, az=az, frange=frange, bins=bins, flags=flags,
+                            apodize_edges=apodize_edges, apodize_edges_samps=apodize_edges_samps,
+                            apodize_flags=apodize_flags, apodize_flags_samps=apodize_flags_samps,)
+    bin_centers = binning_dict['bin_centers']
+    bin_counts = binning_dict['bin_counts']
+    binned_signal = binning_dict['binned_signal']
+    binned_signal_sigma = binning_dict['binned_signal_sigma']
+    uniform_binned_signal_sigma = np.nanmedian(binned_signal_sigma, axis=-1)
+
+    azss_stats = core.AxisManager(aman.dets)
+    azss_stats.wrap('binned_az', bin_centers, [(0, core.IndexAxis('bin_az_samps', count=bins))])
+    azss_stats.wrap('bin_counts', bin_counts, [(0, 'dets'), (1, 'bin_az_samps')])
+    azss_stats.wrap('binned_signal', binned_signal, [(0, 'dets'), (1, 'bin_az_samps')])
+    azss_stats.wrap('binned_signal_sigma', binned_signal_sigma, [(0, 'dets'), (1, 'bin_az_samps')])
+    azss_stats.wrap('uniform_binned_signal_sigma', uniform_binned_signal_sigma, [(0, 'dets')])
+    azss_stats.wrap('method', method)
+    if not frange is None:
+        azss_stats.wrap('frange_min', frange[0])
+        azss_stats.wrap('frange_max', frange[1])
+    if max_mode:
+        azss_stats.wrap('max_mode', max_mode)
+    return azss_stats
+
+def get_azss(aman, signal='signal', az=None, frange=None, bins=100, flags=None,
+            apodize_edges=True, apodize_edges_samps=1600, apodize_flags=True, apodize_flags_samps=200,
             apply_prefilt=True, prefilt_cfg=None, prefilt_detrend='linear',
             method='interpolate', max_mode=None, subtract_in_place=False,
-            merge_stats=True, azss_stats_name='azss_stats',
-            merge_model=True, azss_model_name='azss_model'):
+            merge_stats=True, azss_stats_name='azss_stats', turnaround_info=None,
+            merge_model=True, azss_model_name='azss_model', left_right=False):
     """
     Derive azss (Azimuth Synchronous Signal) statistics and model from the given axismanager data.
     **NOTE:** This function does not modify the ``signal`` unless ``subtract_in_place = True``.
@@ -160,14 +191,14 @@ def get_azss(aman, signal='signal', az=None, range=None, bins=100, flags=None,
         A numpy array representing the signal to be used for azss extraction. If not provided, the signal is taken from aman.signal.
     az: array-like, optional
         A 1D numpy array representing the azimuth angles. If not provided, the azimuth angles are taken from aman.boresight.az.
-    range: list, optional
+    frange: list, optional
         A list specifying the range of azimuth angles to consider for binning. Defaults to [-np.pi, np.pi].
         If None, [min(az), max(az)] will be used for binning.
     bins: int or sequence of scalars
         If bins is an int, it defines the number of equal-width bins in the given range (100, by default).
         If bins is a sequence, it defines the bin edges, including the rightmost edge, allowing for non-uniform bin widths.
-        If `bins` is a sequence, `bins` overwrite `range`.
-    flags : RangesMatrix, optinal
+        If `bins` is a sequence, `bins` overwrite `frange`.
+    flags : RangesMatrix, optional
         Flag indicating whether to exclude flagged samples when binning the signal.
         Default is no mask applied.
     apodize_edges : bool, optional
@@ -189,7 +220,7 @@ def get_azss(aman, signal='signal', az=None, range=None, bins=100, flags=None,
         In 'interpolate', binned signal is used directly.
         In 'fit', fitting is applied to the binned signal.
         Defaults to 'interpolate'.
-    max_mode: integer, optinal
+    max_mode: integer, optional
         The number of Legendre modes to use for azss when method is 'fit'. Required when method is 'fit'.
     subtract_in_place: bool
         If True, it subtract the modeled tod from original signal. The aman.signal will be modified.
@@ -201,6 +232,11 @@ def get_azss(aman, signal='signal', az=None, range=None, bins=100, flags=None,
         Boolean flag indicating whether to merge the azss model with the aman. Defaults to True.
     azss_model_name: string, optional
         The name to assign to the merged azss model. Defaults to 'azss_model'.
+    left_right: bool
+        Default False. If True estimate (and subtract) the AzSS template for left and right subscans
+        separately.
+    turnaround_info: FlagManager or AxisManager
+        Optional, default is aman.flags.
 
     Returns
     -------
@@ -235,43 +271,142 @@ def get_azss(aman, signal='signal', az=None, range=None, bins=100, flags=None,
 
     if az is None:
         az = aman.boresight.az
+
+    if flags is None:
+        flags = Ranges.from_mask(np.zeros(aman.samps.count).astype(bool))
+    else: 
+        flags = aman.flags.reduce(flags=flags, method='union', wrap=False)
+
+    if left_right:
+        if turnaround_info is None:
+            turnaround_info = aman.flags
+        if isinstance(turnaround_info, str):
+            _f = attrgetter(turnaround_info)
+            turnaround_info = _f(aman)
+        if not isinstance(turnaround_info, (core.AxisManager, core.FlagManager)):
+            raise TypeError('turnaround_info must be AxisManager or FlagManager')
+
+        if "valid_right_scans" not in turnaround_info: 
+            left_mask = turnaround_info.left_scan
+            right_mask = turnaround_info.right_scan
+        else:
+            left_mask = turnaround_info.valid_left_scans
+            right_mask = turnaround_info.valid_right_scans
+
+        azss_left = _prepare_azss_stats(aman, signal, az, frange, bins, flags+~left_mask, apodize_edges,
+                                        apodize_edges_samps, apodize_flags, apodize_flags_samps,
+                                        method=method, max_mode=max_mode)
+        azss_left.add_axis(aman.samps)
+        azss_left.add_axis(aman.dets)
+        azss_left.wrap('mask', left_mask, [(0, 'dets'), (1, 'samps')])
+        azss_right = _prepare_azss_stats(aman, signal, az, frange, bins, flags+~right_mask, apodize_edges,
+                                        apodize_edges_samps, apodize_flags, apodize_flags_samps,
+                                        method=method, max_mode=max_mode)
+        azss_right.add_axis(aman.samps)
+        azss_right.add_axis(aman.dets)
+        azss_right.wrap('mask', right_mask, [(0, 'dets'), (1, 'samps')])
         
-    # do binning
-    binning_dict = bin_by_az(aman, signal=signal, az=az, range=range, bins=bins, flags=flags,
-                            apodize_edges=apodize_edges, apodize_edges_samps=apodize_edges_samps, 
-                            apodize_flags=apodize_flags, apodize_flags_samps=apodize_flags_samps,)
-    bin_centers = binning_dict['bin_centers']
-    bin_counts = binning_dict['bin_counts']
-    binned_signal = binning_dict['binned_signal']
-    binned_signal_sigma = binning_dict['binned_signal_sigma']
-    uniform_binned_signal_sigma = np.nanmedian(binned_signal_sigma, axis=-1)
-    
-    azss_stats = core.AxisManager(aman.dets)
-    azss_stats.wrap('binned_az', bin_centers, [(0, core.IndexAxis('bin_az_samps', count=bins))])
-    azss_stats.wrap('bin_counts', bin_counts, [(0, 'dets'), (1, 'bin_az_samps')])
-    azss_stats.wrap('binned_signal', binned_signal, [(0, 'dets'), (1, 'bin_az_samps')])
-    azss_stats.wrap('binned_signal_sigma', binned_signal_sigma, [(0, 'dets'), (1, 'bin_az_samps')])
-    azss_stats.wrap('uniform_binned_signal_sigma', uniform_binned_signal_sigma, [(0, 'dets')])
-    
-    if method == 'fit':
-        if type(max_mode) is not int:
-            raise ValueError('max_mode is not provided as integer')
-        azss_stats, model_sig_tod = fit_azss(az=az, azss_stats=azss_stats, max_mode=max_mode, fit_range=range)
-        
-    if method == 'interpolate':
-        f_template = interp1d(bin_centers, binned_signal, fill_value='extrapolate')
-        model_sig_tod = f_template(aman.boresight.az)
-    
+        azss_stats = core.AxisManager(aman.dets)
+        azss_stats.wrap('azss_stats_left', azss_left)
+        azss_stats.wrap('azss_stats_right', azss_right)
+        azss_stats.wrap('left_right', left_right)
+    else:
+        azss_stats = _prepare_azss_stats(aman, signal, az, frange, bins, flags, apodize_edges,
+                                        apodize_edges_samps, apodize_flags, apodize_flags_samps,
+                                        method=method, max_mode=max_mode)
+        azss_stats.wrap('left_right', left_right)
+
+    model_left, model_right, model = None, None, None
+    if merge_model or subtract_in_place:
+        if left_right:
+            azss_stats, model_left, model_right = get_model_sig_tod(aman, azss_stats, az)
+            if merge_model:
+                aman.wrap(azss_model_name+'_left', model_left, [(0, 'dets'), (1, 'samps')])
+                aman.wrap(azss_model_name+'_right', model_right, [(0, 'dets'), (1, 'samps')])
+            if subtract_in_place:
+                if signal_name is None:
+                    lmask = left_mask.mask()
+                    signal[lmask] -= model_left[lmask].astype(signal.dtype)
+                    rmask = right_mask.mask()
+                    signal[rmask] -= model_right[rmask].astype(signal.dtype)
+                else:
+                    lmask = left_mask.mask()
+                    aman[signal_name][lmask] -= model_left[lmask].astype(aman[signal_name].dtype)
+                    rmask = right_mask.mask()
+                    aman[signal_name][rmask] -= model_right[rmask].astype(aman[signal_name].dtype)
+        else:
+            azss_stats, model, _ = get_model_sig_tod(aman, azss_stats, az)
+            if merge_model:
+                aman.wrap(azss_model_name, model, [(0, 'dets'), (1, 'samps')])
+            if subtract_in_place:
+                if signal_name is None:
+                    signal -= model.astype(signal.dtype)
+                else:
+                    aman[signal_name] -= model.astype(aman[signal_name].dtype)
+
     if merge_stats:
         aman.wrap(azss_stats_name, azss_stats)
-    if merge_model:
-        aman.wrap(azss_model_name, model_sig_tod, [(0, 'dets'), (1, 'samps')])
-    if subtract_in_place:
-        aman[signal_name] = np.subtract(signal, model_sig_tod, dtype='float32')
-    return azss_stats, model_sig_tod
 
-def subtract_azss(aman, signal='signal', azss_template_name='azss_model',
-                  subtract_name='azss_remove', in_place=False, remove_template=True):
+    if left_right:
+        return azss_stats, model_left, model_right
+    else:
+        return azss_stats, model
+
+def get_model_sig_tod(aman, azss_stats, az=None):
+    """
+    Function to return the azss template for subtraction given the azss_stats AxisManager
+    """
+    if az is None:
+        az = aman.boresight.az
+
+    if azss_stats.left_right:
+        model = []
+        for fld in ['azss_stats_left', 'azss_stats_right']:
+            _azss_stats = azss_stats[fld]
+            if 'frange_min' in _azss_stats:
+                frange = (_azss_stats.frange_min, _azss_stats.frange_max)
+            else:
+                frange = None
+            if _azss_stats.method == 'fit':
+                if type(_azss_stats.max_mode) is not int:
+                    raise ValueError('max_mode is not provided as integer')
+
+                _azss_stats, _model = fit_azss(az=az, azss_stats=_azss_stats,
+                                                   max_mode=_azss_stats.max_mode,
+                                                   fit_range=frange)
+                azss_stats.wrap(fld, _azss_stats, overwrite=True)
+                model.append(_model)
+
+            if _azss_stats.method == 'interpolate':
+                good_az = np.logical_and(_azss_stats.binned_az >= np.min(az), _azss_stats.binned_az <= np.max(az))
+                f_template = interp1d(_azss_stats.binned_az[good_az],
+                                      _azss_stats.binned_signal[:, good_az], fill_value='extrapolate')
+                _model = f_template(az)
+                model.append(_model)
+        for ii in range(len(model)):
+            model[ii][~np.isfinite(model[ii])] = 0
+        return azss_stats, model[0], model[1]
+
+    else:
+        if azss_stats.method == 'fit':
+            if type(azss_stats.max_mode) is not int:
+                    raise ValueError('max_mode is not provided as integer')
+            if 'frange_min' in azss_stats:
+                frange = (azss_stats.frange_min, azss_stats.frange_max)
+            else:
+                frange = None
+            azss_stats, model = fit_azss(az=az, azss_stats=azss_stats,
+                                         max_mode=azss_stats.max_mode,
+                                         fit_range=frange)
+        if azss_stats.method == 'interpolate':
+            good_az = np.logical_and(azss_stats.binned_az >= np.min(az), azss_stats.binned_az <= np.max(az))
+            f_template = interp1d(azss_stats.binned_az[good_az], azss_stats.binned_signal[:, good_az], fill_value='extrapolate')
+            model = f_template(az)
+            model[~np.isfinite(model)] = 0
+        return azss_stats, model, None
+
+def subtract_azss(aman, azss_stats, signal='signal', subtract_name='azss_remove',
+                  in_place=False):
     """
     Subtract the scan synchronous signal (azss) template from the
     signal in the given axis manager.
@@ -280,21 +415,17 @@ def subtract_azss(aman, signal='signal', azss_template_name='azss_model',
     ----------
     aman : AxisManager
         The axis manager containing the signal and the azss template.
+    azss_stats: AxisManager
+        Contains AxisManager from get_azss.
     signal : str, optional
         The name of the field in the axis manager containing the signal to be processed.
         Defaults to 'signal'.
-    azss_template_name : str, optional
-        The name of the field in the axis manager containing the azss template.
-        Defaults to 'azss_model'.
     subtract_name : str, optional
         The name of the field in the axis manager that will store the azss-subtracted signal.
         Only used if in_place is False. Defaults to 'azss_remove'.
     in_place : bool, optional
         If True, the subtraction is done in place, modifying the original signal in the axis manager.
         If False, the result is stored in a new field specified by subtract_name. Defaults to False.
-    remove_template : bool, optional
-        If True, the azss template field is removed from the axis manager after subtraction.
-        Defaults to True.
 
     Returns
     -------
@@ -313,15 +444,31 @@ def subtract_azss(aman, signal='signal', azss_template_name='azss_model',
     else:
         raise TypeError("Signal must be None, str, or ndarray")
 
+    azss_stats, model_left, model_right = get_model_sig_tod(aman, azss_stats, az=None)
+
     if in_place:
         if signal_name is None:
-            signal -= aman[azss_template_name].astype(signal.dtype)
+            if azss_stats.left_right:
+                for model, azss_fld in zip([model_left, model_right], ['azss_stats_left', 'azss_stats_right']):
+                    mask = azss_stats[azss_fld]['mask'].mask()
+                    signal[:, mask] -= model[:, mask].astype(signal.dtype)
+            else:
+                signal -= model_left.astype(signal.dtype)
         else:
-            aman[signal_name] -= aman[azss_template_name].astype(aman[signal_name].dtype)
+            if azss_stats.left_right:
+                for model, azss_fld in zip([model_left, model_right], ['azss_stats_left', 'azss_stats_right']):
+                    mask = azss_stats[azss_fld]['mask'].mask()
+                    aman[signal_name][mask] -= model[mask].astype(aman[signal_name].dtype)
+            else:
+                aman[signal_name] -= model_left.astype(aman[signal_name].dtype)
     else:
-        aman.wrap(subtract_name, 
-                  np.subtract(aman[signal_name], aman[azss_template_name], dtype='float32'),
-                  [(0, 'dets'), (1, 'samps')])
-        
-    if remove_template:
-        aman.move(azss_template_name, None)
+        if azss_stats.left_right:
+            wrap_sig = np.copy(signal)
+            for model, azss_fld in zip([model_left, model_right], ['azss_stats_left', 'azss_stats_right']):
+                mask = azss_stats[azss_fld]['mask'].mask()
+                wrap_sig[:, mask] -= model[:, mask].astype(signal.dtype)
+            aman.wrap(subtract_name, wrap_sig, [(0, 'dets'), (1, 'samps')])
+        else:
+            aman.wrap(subtract_name, 
+                      np.subtract(aman[signal_name], model_left, dtype='float32'),
+                      [(0, 'dets'), (1, 'samps')])
diff --git a/sotodlib/tod_ops/deproject.py b/sotodlib/tod_ops/deproject.py
index feb3271ca..c05847b8a 100644
--- a/sotodlib/tod_ops/deproject.py
+++ b/sotodlib/tod_ops/deproject.py
@@ -49,7 +49,7 @@ def get_qu_common_mode_coeffs(aman, Q_signal=None, U_signal=None, merge=False):
     return output_aman
 
 def subtract_qu_common_mode(aman, Q_signal=None, U_signal=None, coeff_aman=None,
-                            merge=False):
+                            merge=False, subtract=True):
     """
     Subtracts the median signal (template) from each detector scaled by the a 
     coupling coefficient per detector.
diff --git a/sotodlib/tod_ops/flags.py b/sotodlib/tod_ops/flags.py
index ed7392927..daf71290c 100644
--- a/sotodlib/tod_ops/flags.py
+++ b/sotodlib/tod_ops/flags.py
@@ -14,7 +14,6 @@
 from .. import coords
 from . import filters
 from . import fourier_filter 
-from . import sub_polyf
 
 def get_det_bias_flags(aman, detcal=None, rfrac_range=(0.1, 0.7),
                        psat_range=None, rn_range=None, si_nan=False,
@@ -718,7 +717,7 @@ def get_inv_var_flags(aman, signal_name='signal', nsigma=5,
 
     return mskinvar
 
-def get_subscans(aman, merge=True, include_turnarounds=False, overwrite=True):
+def get_subscans(aman, flags="flags", merge=True, include_turnarounds=False, overwrite=True):
     """
     Returns an axis manager with information about subscans.
     This includes direction and a ranges matrix (subscans samps)
@@ -744,6 +743,7 @@ def get_subscans(aman, merge=True, include_turnarounds=False, overwrite=True):
     """
     if not include_turnarounds:
         ss_ind = (~aman.flags.turnarounds).ranges() # sliceable indices (first inclusive, last exclusive) for subscans
+        #ss_ind = (~aman[flags]["turnarounds"]).ranges()
     else:
         left = aman.flags.left_scan.ranges()
         right = aman.flags.right_scan.ranges()
@@ -998,3 +998,81 @@ def noise_fit_flags(aman, low_wn, high_wn, high_fk):
         return flag_valid_fk
     else:
         return None
+
+def get_noisy_subscan_flags(aman, subscan_stats, flags="flags",
+                         nstd_lim=None, ptp_lim=None, kurt_lim=None, 
+                         skew_lim=None, noisy_detector_lim=0.9,
+                         merge=False, overwrite=False, 
+                         name="noisy_subscan"):
+    """
+    Identify and flag bad subscans based on various statistical thresholds.
+    aman : AxisManager 
+        The tod.
+    subscan_stats : dict
+        Dictionary containing statistical metrics for subscans. Keys should 
+        include 'std', 'ptp', 'kurtosis', and 'skew'.
+    nstd_lim : float [optional]
+        Threshold for standard deviation.
+    ptp_lim : float [optional]
+        Threshold for peak-to-peak values in pW.
+    kurt_lim : float [optional]
+        Threshold for kurtosis.
+    skew_lim : float [optional]
+        Threshold for skewness.
+    merge : bool
+        If true, merges the generated flag into aman.
+    overwrite : bool
+        If true, write over flag. If false, don't.
+    name : str
+        Name of flag to add to aman.flags if merge is True.
+
+    Returns
+    -------
+    noisy_subscan_flags : RangesMatrix
+        RangesMatrix of bad subscans.
+    noisy_detector_flags : ndarray
+        Array indicating detectors that are too noisy for more than 50% of the subscan duration.
+
+    """
+
+    if 'flags' not in aman:
+        overwrite = False
+        merge = False
+    if overwrite and name in aman.flags:
+        aman.flags.move(name, None)
+    
+    median_std = np.median(subscan_stats["std"], axis=1)[:, np.newaxis]
+    noisy_subscan_indicator = np.zeros_like(subscan_stats["std"], dtype=bool)
+    
+    if ptp_lim is not None and 'ptp' in subscan_stats:
+        noisy_subscan_indicator |= (subscan_stats['ptp'] > ptp_lim)
+    
+    if nstd_lim is not None and 'std' in subscan_stats:
+        noisy_subscan_indicator |= (subscan_stats['std'] > median_std * nstd_lim)
+    
+    if kurt_lim is not None and 'kurtosis' in subscan_stats:
+        noisy_subscan_indicator |= (np.abs(subscan_stats['kurtosis']) > kurt_lim)
+    
+    if skew_lim is not None and 'skew' in subscan_stats:
+        noisy_subscan_indicator |= (np.abs(subscan_stats['skew']) > skew_lim)
+
+    if "subscan_info" in aman:
+        ssf = aman.subscan_info.subscan_flags
+    elif "subscans" in aman: 
+        # If reading from proc_aman
+        ssf = aman.subscans
+    else:
+        get_subscans(aman, flags=flags, merge=False)
+
+    zeros=RangesMatrix.zeros((1, aman.samps.count)) # Needed when no subscans are True
+    # For each det, select flagged subscans and merge their Ranges
+    noisy_subscan_flags = RangesMatrix([np.sum(RangesMatrix.concatenate([ssf[noisy_subscan_indicator[idet]], zeros]), axis=0) for idet in range(aman.dets.count)])
+    # Detectors which are too noisy for > 50% of the subscan duration
+    noisy_detector_flags = np.array([np.sum(np.diff(rng.ranges(), axis=1)) / aman.samps.count for rng in noisy_subscan_flags.ranges]) > noisy_detector_lim
+    
+    if merge:
+        if name in aman.flags and not overwrite:
+            raise ValueError(f"Flag name {name} already exists in aman.flags")
+        aman.flags.wrap(name, noisy_subscan_flags)
+
+    return noisy_subscan_flags, ~noisy_detector_flags
diff --git a/tests/test_azss.py b/tests/test_azss.py
index e3eb779a2..66fa0c8af 100644
--- a/tests/test_azss.py
+++ b/tests/test_azss.py
@@ -104,7 +104,7 @@ class AzssTest(unittest.TestCase):
     def test_fit(self):
         max_mode = 10
         tod = make_fake_azss_tod(noise_amp=0, n_scans=50, max_mode=max_mode)
-        azss_stats, model_sig_tod = azss.get_azss(tod, method='fit', max_mode=max_mode, range=None, bins=100)
+        azss_stats, model_sig_tod = azss.get_azss(tod, method='fit', max_mode=max_mode, frange=None, bins=100)
         ommedian = get_coeff_metric(tod)
         print(ommedian)
         self.assertTrue(ommedian < 5.0)