Repository: spacetelescope/jwst

Issue:
Implementation of 1/f correction
_Issue [JP-3639](https://jira.stsci.edu/browse/JP-3639) was created on JIRA by [David Law](https://jira.stsci.edu/secure/ViewProfile.jspa?name=dlaw):_

Various other tickets exist describing the need for 1/f correction, e.g., <https://jira.stsci.edu/browse/JP-2576> discussing the rationale for doing so.  Likewise, <https://outerspace.stsci.edu/pages/viewpage.action?spaceKey=JWSTCC&title=JWST+Pipeline+1-over-f+Noise+Removal> contains notes from various meetings discussing potential algorithms.

This ticket focuses on establishing the detailed method of implementation in the pipeline, and is separate to more clearly distinguish discussion such details from discussion of the core algorithms wrapped by the step.

The new '1/f correction' step in the pipeline should be immediately prior to the current ramp_fitting step in calwebb_detector1.  The reason for this is that it is best to correct 1/f noise at the group stage in order to reliably correct the 1/f signal both in blank regions and underneath sources.

This 1/f correction would be done using the source-free science pixels in the detector using some TBD algorithm.  Regardless of the details of this algorithm though, it will be necessary to use a mask to mask out sources in the scene, measure the 1/f signal, and apply the correction to the group stage data prior to ramp fitting.  This step should be capable of using such a masked scene if passed in to the step, and should otherwise attempt to compute such a mask on the fly.  In order to do so, the 1/f step would need to call ramp_fitting to produce a 'preliminary' rate image with which to do the masking.  As such ramp fitting could strictly be done twice- once within the 1/f step to generate the on-the-fly mask, and then again in the regular ramp_fitting stage using the 1/f corrected ramp.  If doing so, the ramp fitting call within the 1/f step should probably use the C-based implementation (or a comparably fast one) by default.

Note that I've identified all four instruments in this ticket.  For NIRCam and NIRISS, this is because the work is most relevant to them.  NIRSpec already has a solution in NSClean, but the proposal here is to move that step into the infrastructure proposed here to have a single location for 1/f correction in the pipeline.  MIRI does not experience 1/f noise, but nonetheless has vertical striping with some similarities to 1/f noise (resulting from variations in the detector dark) that may be possible to correct in a similar manner and place in the pipeline.

Additional details:
 * This step should be inserted between the 'jump' and 'ramp_fitting' steps.  It should be set to skip by default for the time being (and maybe permanently, TBD).
 * This step should read in a single group-level file and an optional 2d mask with x/y pixel size the same as the input group-data file, and write out modified group-level data.  No reference files for this step at the current time.
 * Within the code it should look to see if an input 2d mask was provided, and if not calculate one from the group data by running it through a temporary ramp fitting (using the C-based implementation) and applying some TBD algorithm to derive the mask.  The 1/f tiger team will fill in details on this masking algorithm in the next 1-2 weeks.
 * Once a mask exists, the group-stage data should be corrected using a choice of two possible algorithms.  The first is the NSClean algorithm, already implemented for rate-stage data, and the second is Chris Willott's 1/f removal algorithm ([https://jwst-docs.stsci.edu/files/255987666/255987671/1/1710512651770/example_nrc_image1overf_5.ipynb]).  Note that this will need to be tweaked slightly to not subtract off the overall increasing baseline counts from the background signal, only the differential counts across the image.  The 1/f team will weigh in on details in the next 1-2 weeks.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
Attention: Patch coverage is `12.53823%` with `286 lines` in your changes missing coverage. Please review.
> Project coverage is 59.99%. Comparing base [(`ccf871a`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/ccf871a479dad57b0a1dc235ea385ed145325d1f?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) to head [(`9a7ad8c`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/9a7ad8c6b7626cd7b9d2899d2b91fdf40a756b4a?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).
> Report is 20 commits behind head on master.

| [Files](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Patch % | Lines |
|---|---|---|
| [jwst/clean\_noise/clean\_noise.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?src=pr&el=tree&filepath=jwst%2Fclean_noise%2Fclean_noise.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9jbGVhbl9ub2lzZS9jbGVhbl9ub2lzZS5weQ==) | 8.78% | [270 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |
| [jwst/clean\_noise/clean\_noise\_step.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?src=pr&el=tree&filepath=jwst%2Fclean_noise%2Fclean_noise_step.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9jbGVhbl9ub2lzZS9jbGVhbl9ub2lzZV9zdGVwLnB5) | 44.44% | [10 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |
| [jwst/nsclean/nsclean\_step.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?src=pr&el=tree&filepath=jwst%2Fnsclean%2Fnsclean_step.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9uc2NsZWFuL25zY2xlYW5fc3RlcC5weQ==) | 16.66% | [5 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |
| [jwst/pipeline/calwebb\_detector1.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?src=pr&el=tree&filepath=jwst%2Fpipeline%2Fcalwebb_detector1.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9waXBlbGluZS9jYWx3ZWJiX2RldGVjdG9yMS5weQ==) | 50.00% | [1 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #8669      +/-   ##
==========================================
+ Coverage   59.56%   59.99%   +0.42%     
==========================================
  Files         391      372      -19     
  Lines       39285    38807     -478     
==========================================
- Hits        23402    23283     -119     
+ Misses      15883    15524     -359     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/spacetelescope/jwst/pull/8669?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

melanieclarke: Regression tests started here:
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1707/
melanieclarke: In the regression test results:
- Minor changes are expected for nsclean for NIRSpec spec2 and downstream products.  I've inspected the differences locally and they are small, mostly due to some revisions in how the scene mask is created, including a new exclusion for large negative outliers.
- Changes from truth for the new regtest for MIRI clean_flicker_noise were just because of differences between my system and Jenkins: my local CRDS cache had the wrong emicorr reference file.  I okified those changes.
- Other changes (NIRISS rate and MIRI badpix_selfcal) are unrelated.

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index b9f075c9c..51fdfe145 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -42,6 +42,18 @@ badpix_selfcal
 - Subtract pedestal dark when constructing min array across selfcal exposures
   for MIRI MRS data. [#8771]
 
+calwebb_detector1
+-----------------
+
+- Added the optional ``clean_flicker_noise`` step between ``jump`` and
+  ``ramp_fit``. [#8669]
+
+clean_flicker_noise
+-------------------
+
+- Implemented this new optional step to clean transient flicker noise (e.g. 1/f noise)
+  from group images in ramp data. [#8669]
+
 cube_build
 ----------
 
@@ -115,6 +127,10 @@ nsclean
   Due to the computational costs of matrix operations, this is a large
   speedup that has little effect on the results. [#8745]
 
+- Merged implementation with the new ``clean_flicker_noise`` step. This step
+  can still be called from the ``calwebb_spec2`` pipeline on NIRSpec rate
+  data, but it is now deprecated. [#8669]
+
 outlier_detection
 -----------------
 
diff --git a/docs/jwst/clean_flicker_noise/arguments.rst b/docs/jwst/clean_flicker_noise/arguments.rst
new file mode 100644
index 000000000..8fc131709
--- /dev/null
+++ b/docs/jwst/clean_flicker_noise/arguments.rst
@@ -0,0 +1,71 @@
+.. _clean_flicker_noise_arguments:
+
+Step Arguments
+==============
+
+The ``clean_flicker_noise`` step has the following optional arguments to control
+the behavior of the processing.
+
+``--fit_method`` (str, default='median')
+  The noise fitting algorithm to use.  Options are 'fft' and 'median'.
+
+``--fit_by_channel`` (boolean, default=False)
+  If set, flicker noise is fit independently for each detector channel.
+  Ignored for MIRI, for subarray data, and for `fit_method` = 'fft'.
+
+``--background_method`` (str, default='median')
+  If 'median', the preliminary background to remove and restore
+  after fitting the noise is a simple median of the background data.
+  If 'model', the background data is fit with a low-resolution model
+  via `~photutils.background.Background2D`.
+  If None, the background value is set to 0.0.
+
+``--background_box_size`` (list of int, default=(32,32))
+  Box size for the data grid used by `~photutils.background.Background2D`
+  when `background_method` = 'model'. For best results, use a
+  box size that evenly divides the input image shape.
+
+``--mask_science_regions`` (boolean, default=False)
+  For NIRSpec, mask regions of the image defined by WCS bounding
+  boxes for slits/slices, as well as any regions known to be
+  affected by failed-open MSA shutters.  For MIRI imaging, mask
+  regions of the detector not used for science.
+
+``--n_sigma`` (float, default=2.0)
+  The sigma-clipping threshold to use when searching for outliers
+  and illuminated pixels to be excluded from use in the background
+  and noise fitting processes.
+
+``--fit_histogram`` (boolean, default=False)
+  If set, the 'sigma' used with `n_sigma` for clipping outliers
+  is derived from a Gaussian fit to a histogram of values.
+  Otherwise, a simple iterative sigma clipping is performed.
+
+``--single_mask`` (boolean, default=True)
+  If set, a single mask will be created, regardless of
+  the number of input integrations. Otherwise, the mask will
+  be a 3D cube, with one plane for each integration.
+
+``--user_mask`` (string, default=None)
+  Path to a user-supplied mask file. If supplied, the mask is used
+  directly and the process of creating a scene mask in the step is
+  skipped.
+
+  The mask file must contain either a `~jwst.datamodels.ImageModel`
+  or a `~jwst.datamodels.CubeModel`, with image dimensions matching
+  the input science data.  If an ImageModel is provided, the same
+  mask will be used for all integrations.  If a CubeModel is provided,
+  the number of slices must equal the number of integrations in
+  the input science data.
+
+``--save_mask`` (boolean, default=False)
+  If set, the mask constructed by the step will be saved to a file
+  with suffix 'mask'.
+
+``--save_background`` (boolean, default=False)
+  If set, the background fit to the group diff images will be saved
+  to a file with suffix 'flicker_bkg'.
+
+``--save_noise`` (boolean, default=False)
+  If set, the residual noise fit and removed from the input data
+  will be saved to a file with suffix 'flicker_noise'.
diff --git a/docs/jwst/clean_flicker_noise/index.rst b/docs/jwst/clean_flicker_noise/index.rst
new file mode 100644
index 000000000..aa8b5e426
--- /dev/null
+++ b/docs/jwst/clean_flicker_noise/index.rst
@@ -0,0 +1,13 @@
+.. _clean_flicker_noise_step:
+
+===================
+Clean Flicker Noise
+===================
+.. toctree::
+   :maxdepth: 2
+
+   main.rst
+   reference_files.rst
+   arguments.rst
+   
+.. automodapi:: jwst.clean_flicker_noise
diff --git a/docs/jwst/clean_flicker_noise/main.rst b/docs/jwst/clean_flicker_noise/main.rst
new file mode 100644
index 000000000..c7356c690
--- /dev/null
+++ b/docs/jwst/clean_flicker_noise/main.rst
@@ -0,0 +1,253 @@
+Description
+===========
+
+:Class: `jwst.clean_flicker_noise.CleanFlickerNoiseStep`
+:Alias: clean_flicker_noise
+
+Overview
+--------
+The ``clean_flicker_noise`` step removes flicker noise from calibrated
+ramp images, after the :ref:`jump <jump_step>` step and prior to
+performing the :ref:`ramp_fitting <ramp_fitting_step>` step.
+
+For NIR detectors, the noise addressed by this step is 1/f noise, which
+appears as faint banding along the detector slow axis.  For NIRSpec and
+NIRISS, 1/f noise looks like vertical stripes; for NIRCam, it appears
+as horizontal stripes.
+
+MIRI images also often show a faint vertical banding, similar to 1/f noise
+but from a different physical source.  This type of flicker noise can be
+corrected by similar methods, so it is also addressed by this step.
+
+To correct for flicker noise, the algorithm requires that the noise
+generated between one group readout and the next be isolated as much
+as possible from the astronomical signal. The residual noise may then
+be fit in frequency space, via an FFT, or may be characterized by a
+median along the rows or columns, as appropriate for the detector.
+The modeled noise is then directly subtracted from each group readout.
+
+The correction is available for any type of NIRSpec, NIRISS, or NIRCam
+exposure. For MIRI, it is available only for imaging exposures.
+
+Creation of a scene mask
+------------------------
+One of the key components of the correction is knowing which pixels can
+be used to fit the background noise.  The step builds a scene mask
+on the fly from a draft rate image, generated from the input ramp data,
+which is used to mark usable and unusable pixels. The mask is a 2D
+Boolean array, having the same size as the image, with
+pixels set to True interpreted as being OK to use.
+
+The process of building the mask varies somewhat depending on the
+observing mode of the image being processed. Some features are common
+to all modes, while others are mode-specific. The following sections
+describe each type of masking that can be applied. At the end, there
+is a summary of the types applied to each instrument mode.
+
+The user-settable step parameter `save_mask` can be used to save the
+scene mask to a file, if desired (see the
+:ref:`step arguments <clean_flicker_noise_arguments>`).
+
+Note that a user may also supply their own mask image as input to the step,
+in which case the process of creating a mask is skipped. The step parameter
+`user_mask` is used to specify an input mask.  If specified, the input
+mask must contain a datamodel matching the shape of a 'rate' or 'rateints'
+file generated for the input science data (either `~jwst.datamodels.ImageModel`
+or `~jwst.datamodels.CubeModel`).  To generate a custom mask, it may be
+easiest to save the mask output from a default run of the ``clean_flicker_noise``
+step on the input data, then modify it as needed.
+
+NIRSpec IFU Slices
+^^^^^^^^^^^^^^^^^^
+For IFU images, the majority of the mask is based on knowing which
+pixels are contained within the footprints of the IFU slices. To do
+this, the image's World Coordinate System (WCS) object is queried in
+order to determine which pixels are contained within each of the 30
+slices. Pixels within each slice are set to False (do not use) in the
+mask.
+
+NIRSpec MOS/FS Slits
+^^^^^^^^^^^^^^^^^^^^
+The footprints of each open MOS slitlet or fixed slit are flagged in
+a similar way as IFU slices. For MOS and FS images, the WCS object is
+queried to determine which pixels are contained within each open
+slit/slitlet and they are set to False in the mask.
+
+NIRSpec MSA Failed Open Shutters
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Pixels affected by stuck open MSA shutters are masked, because they
+may contain signal. This is accomplished by setting all pixels flagged by the
+:ref:`msaflagopen <msaflagopen_step>` step with DQ value "MSA_FAILED_OPEN"
+to False in the mask.
+
+NIRSpec Fixed-Slit Region Pixels
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Full-frame MOS and IFU images may contain signal from the always open
+fixed slits, which appear in a fixed region in the middle of each image.
+The entire region containing the fixed slits is masked out when
+processing MOS and IFU images. The masked region is currently hardwired
+in the step to image indexes [1:2048, 923:1116], where the indexes are
+in x, y order and in 1-indexed values.
+
+Note, however, that it is possible to plan one or more fixed slit targets
+alongside MSA slitlets in MOS observations. In this situation, the fixed
+slit region is not automatically masked.
+
+MIRI Imaging Metering Structure
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+The MIRI imager has a metering structure covering a large part of the
+detector area. These regions must be masked in order to fit and
+remove the background data in the science areas of the detector.
+Regions marked with DQ value "DO_NOT_USE" by the
+:ref:`flat_field <flatfield_step>` step are set to False in the
+scene mask.
+
+Missing Data
+^^^^^^^^^^^^
+Any pixel in the draft rate image that has a value of NaN or exactly zero
+is flagged as False in the mask. This typically includes any reference
+pixels that are present in the exposure.
+
+Outliers
+^^^^^^^^
+For imaging modes, bright, compact sources must be distinguished
+from the background and masked in order to fit and remove the
+smooth background level.
+
+For spectral modes that already have significant masking applied,
+pixels in the unilluminated areas of the region can still contain anomalous
+signal, due to uncaught cosmic rays, hot pixels, etc.
+
+For both modes, a sigma-clipping routine is employed to find such outliers
+within the draft rate image and set them to False in the mask. All pixels with
+values greater than :math:`median+n\_sigma*sigma` are assumed to contain
+signal and are set to False in the scene mask. In addition, all pixels
+with values less than :math:`median-3.0*sigma` are assumed to be bad pixels,
+and are also set to False in the scene mask.
+
+Mode-Specific Masking Steps
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+The following table indicates which flavors of masking are applied to
+images from each instrument and observing mode.
+
+.. |c| unicode:: U+2713 .. checkmark
+
++--------------------------+-----+-----+-----+-------+--------+--------+
+|                          |    NIRSpec      | MIRI  | NIRCam | NIRISS |
++--------------------------+-----+-----+-----+-------+--------+--------+
+|                          | IFU | MOS |  FS | Image | All    | All    |
++==========================+=====+=====+=====+=======+========+========+
+| IFU Slices\ :sup:`1`     | |c| |     |     |       |        |        |
++--------------------------+-----+-----+-----+-------+--------+--------+
+| Slits/Slitlets\ :sup:`1` |     | |c| | |c| |       |        |        |
++--------------------------+-----+-----+-----+-------+--------+--------+
+| MSA_FAILED_OPEN\ :sup:`1`| |c| | |c| |     |       |        |        |
++--------------------------+-----+-----+-----+-------+--------+--------+
+| Non-science\ :sup:`1`    |     |     |     | |c|   |        |        |
++--------------------------+-----+-----+-----+-------+--------+--------+
+| FS Region\ :sup:`1`      | |c| | |c| |     |       |        |        |
++--------------------------+-----+-----+-----+-------+--------+--------+
+| Missing Data             | |c| | |c| | |c| | |c|   | |c|    | |c|    |
++--------------------------+-----+-----+-----+-------+--------+--------+
+| Outliers                 | |c| | |c| | |c| | |c|   | |c|    | |c|    |
++--------------------------+-----+-----+-----+-------+--------+--------+
+
+:sup:`1`\ These steps are only applied if the
+:ref:`step parameter <clean_flicker_noise_arguments>`
+`mask_science_regions` is set to True.
+
+Correction Algorithm
+--------------------
+
+The detailed process for fitting and removing flicker noise is as follows.
+See the :ref:`step arguments <clean_flicker_noise_arguments>` for more
+information on all referenced parameters.
+
+#. From the calibrated ramp input, make a draft rate (`single_mask` = True)
+   or rateints (`single_mask` = False) file.
+
+#. Create a scene mask from the rate data.
+
+   #. If `mask_science_regions` is set and the input is NIRSpec data,
+      run :ref:`assign_wcs <assign_wcs_step>` and
+      :ref:`msaflagopen <msaflagopen_step>` on the draft rate data,
+      then mask any known science areas or failed-open MSA shutters.
+
+      This will mask out regions that are likely to contain significant
+      astronomical signal.
+
+   #. If `mask_science_regions` is set and the input is MIRI imaging data,
+      run :ref:`flat_field <flatfield_step>` on the draft rate data,
+      and extract just the DQ plane from the output. Pixels flagged
+      as 'DO_NOT_USE' by the flat fielding process are masked.
+
+      This will mask out regions of the detector under the metering
+      structure.
+
+   #. Iteratively sigma clip the data to get a center value (mean or median)
+      and sigma value (standard deviation).
+
+   #. If `fit_histogram` is set, compute a histogram from 4-sigma clipped
+      values and fit a Gaussian to it to refine the center and sigma values.
+
+   #. Mask data more than 3 * sigma below the center as bad values.
+
+   #. Mask data more than `n_sigma` * sigma above the center as signal
+      (not background).
+
+#. Iterate over each integration and group in the data, to fit and correct
+   for noise.
+
+   #. Make a diff image (current group – previous group) to correct.
+
+   #. Fit and remove a background level, using the scene mask to identify
+      background pixels.
+
+      #. Clip the background data in the diff image to remove more outliers.
+
+      #. If `background_method` = 'median', the background value is a simple
+         median of the remaining values.
+
+      #. If `background_method` = 'model', the background data is fit with
+         a low-resolution model via the photutils
+         `Background2D <https://photutils.readthedocs.io/en/latest/api/photutils.background.Background2D.html>`_
+         utility. The resolution box size is set by `background_box_size`.
+
+      #. Subtract the background level from the diff image and clip again
+         to `n_sigma` * sigma, with sigma recomputed from the
+         background-subtracted data in the remaining background pixels.
+
+   #. Fit and remove the residual noise in the background-subtracted image.
+
+      #. If `fit_method` = 'fft', the ``nsclean`` library is called to fit
+         and remove the noise in frequency space.
+
+      #. If `fit_method` = 'median', the noise is fit with a simple median
+         along the appropriate detector axis and subtracted from the
+         background-subtracted image.
+
+         If `fit_by_channel` = True, and the data is a NIR full-frame exposure,
+         the median value is computed and subtracted independently for each
+         detector channel.
+
+   #. Restore the background level to the cleaned, background-subtracted
+      diff image.
+
+   #. Add the cleaned diff back to a cleaned version of the previous
+      group image.
+
+References
+==========
+
+The FFT cleaning algorithm implementation is based on NSClean,
+developed by Bernard Rauscher. Details on the source of the correlated
+noise and the algorithm used by the ``nsclean`` library to fit and
+remove it can be found in
+`Rauscher 2024 <https://ui.adsabs.harvard.edu/abs/2023arXiv230603250R/abstract>`__.
+
+The background fitting and median cleaning algorithm are based on
+the image1overf algorithm, developed by Chris Willott, and available
+on GitHub at `chriswillott/jwst <https://github.com/chriswillott/jwst>`__.
+The algorithm was adapted to the `clean_flicker_noise` step and is released
+under the BSD license for the JWST calibration pipeline by permission
+of the author.
diff --git a/docs/jwst/clean_flicker_noise/reference_files.rst b/docs/jwst/clean_flicker_noise/reference_files.rst
new file mode 100644
index 000000000..552e1717a
--- /dev/null
+++ b/docs/jwst/clean_flicker_noise/reference_files.rst
@@ -0,0 +1,3 @@
+Reference Files
+===============
+The ``clean_flicker_noise`` step does not use any reference files.
diff --git a/docs/jwst/nsclean/arguments.rst b/docs/jwst/nsclean/arguments.rst
index c91bfaba3..9772ce649 100644
--- a/docs/jwst/nsclean/arguments.rst
+++ b/docs/jwst/nsclean/arguments.rst
@@ -6,19 +6,65 @@ Step Arguments
 The ``nsclean`` step has the following optional arguments to control
 the behavior of the processing.
 
+``--fit_method`` (str, default='fft')
+  The noise fitting algorithm to use.  Options are 'fft' and 'median'.
+
+``--fit_by_channel`` (boolean, default=False)
+  If set, flicker noise is fit independently for each detector channel.
+  Ignored for subarray data and for `fit_method` = 'fft'.
+
+``--background_method`` (str, default=None)
+  If 'median', the preliminary background to remove and restore
+  after fitting the noise is a simple median of the background data.
+  If 'model', the background data is fit with a low-resolution model
+  via `~photutils.background.Background2D`.
+  If None, the background value is set to 0.0.
+
+``--background_box_size`` (list of int, default=(32,32))
+  Box size for the data grid used by `~photutils.background.Background2D`
+  when `background_method` = 'model'. For best results, use a
+  box size that evenly divides the input image shape.
+
 ``--mask_spectral_regions`` (boolean, default=True)
-  Mask regions in IFU and MOS images that are within the bounding boxes
-  for each slice or slitlet defined in the WCS object of the image.
+  Mask regions of the image defined by WCS bounding
+  boxes for slits/slices, as well as any regions known to be
+  affected by failed-open MSA shutters.
 
 ``--n_sigma`` (float, default=5.0)
   The sigma-clipping threshold to use when searching for outliers
-  and illuminated pixels to be excluded from use in the fitting
-  process.
+  and illuminated pixels to be excluded from use in the background
+  and noise fitting processes.
 
-``--save_mask`` (boolean, default=False)
-  A flag to indicate whether the mask constructed by the step
-  should be saved to a file.
+``--fit_histogram`` (boolean, default=False)
+  If set, the 'sigma' used with `n_sigma` for clipping outliers
+  is derived from a Gaussian fit to a histogram of values.
+  Otherwise, a simple iterative sigma clipping is performed.
+
+``--single_mask`` (boolean, default=False)
+  If set, a single mask will be created, regardless of
+  the number of input integrations. Otherwise, the mask will
+  be a 3D cube, with one plane for each integration.
 
 ``--user_mask`` (string, default=None)
   Path to a user-supplied mask file. If supplied, the mask is used
-  directly and the process of creating a mask in the step is skipped.
+  directly and the process of creating a scene mask in the step is
+  skipped.
+
+  The mask file must contain either a `~jwst.datamodels.ImageModel`
+  or a `~jwst.datamodels.CubeModel`, with image dimensions matching
+  the input science data.  If an ImageModel is provided, the same
+  mask will be used for all integrations.  If a CubeModel is provided,
+  the number of slices must equal the number of integrations in
+  the input science data.
+
+``--save_mask`` (boolean, default=False)
+  If set, the mask constructed by the step will be saved to a file
+  with suffix 'mask'.
+
+``--save_background`` (boolean, default=False)
+  If set, the background fit to the group diff images will be saved
+  to a file with suffix 'flicker_bkg'.
+
+``--save_noise`` (boolean, default=False)
+  If set, the residual noise fit and removed from the input data
+  will be saved to a file with suffix 'flicker_noise'.
diff --git a/docs/jwst/nsclean/main.rst b/docs/jwst/nsclean/main.rst
index 161e84d9d..40dcda7a6 100644
--- a/docs/jwst/nsclean/main.rst
+++ b/docs/jwst/nsclean/main.rst
@@ -4,152 +4,26 @@ Description
 :Class: `jwst.nsclean.NSCleanStep`
 :Alias: nsclean
 
-Overview
-========
-The ``nsclean`` step applies an algorithm for removing correlated read
-noise from NIRSpec images. The noise often appears as faint vertical
-banding and so-called "picture frame noise." The algorithm uses dark
-(unilluminated) areas of an image to fit a background model in Fourier
-space. When the fit is subtracted, it removes nearly all correlated noise.
-Compared to simpler strategies, like subtracting a rolling median, this
-algorithm is more thorough and uniform. It is also computationally
-undemanding, typically requiring only a few seconds to clean a full-frame
-image.
-
-The correction can be applied to any type of NIRSpec exposure, including
-IFU, MOS, fixed slit, and Bright Object Time Series (BOTS), in both full-frame
-and subarray readouts. Time series (3D) data are corrected one integration
-at a time.
-
-Details on the source of the correlated noise and the algorithm used
-in the ``nsclean`` step to fit and remove it can be found in
-`Rauscher 2023 <https://ui.adsabs.harvard.edu/abs/2023arXiv230603250R/abstract>`_.
-
-Upon completion of the step, the step status keyword "S_NSCLEN" gets set
-to "COMPLETE" in the output science data.
-
-Assumptions
-===========
-As described below, the creation of a pixel mask depends on the presence
-of a World Coordinate System (WCS) object for the image, which is
-constructed by the :ref:`assign_wcs <assign_wcs_step>` step.
-In addition, creating a mask for IFU and MOS images depends on
-the presence of DQ flags assigned by the
-:ref:`msaflagopen <msaflagopen_step>` step.
-It is therefore required that those steps be run before attempting to
-apply ``nsclean``.
-
-Creation of an image mask
-=========================
-One of the key components of the correction is knowing which pixels are
-unilluminated and hence can be used in fitting the background noise.
-The step builds a mask on the fly for each image, which is used to mark
-useable and unuseable pixels. The mask is a 2D boolean array, having the same
-size as the image, with pixels set to True interpreted as being OK to use.
-The process of building the mask varies somewhat depending on the
-observing mode of the image being processed. Some features are common
-to all modes, while others are mode-specific. The following sections
-describe each type of masking that can be applied and at the end there
-is a summary of the types applied to each image mode.
-
-The user-settable step parameter `save_mask` can be used to save the
-mask to a file, if desired (see :ref:`nsclean step arguments <nsclean_arguments>`).
-
-Note that a user may supply their own mask image as input to the step,
-in which case the process of creating a mask is skipped. The step parameter
-`user_mask` is used to specify an input mask.
-
-IFU Slices
-----------
-For IFU images the majority of the mask is based on knowing which
-pixels are contained within the footprints of the IFU slices. To do
-this, the image's World Coordinate System (WCS) object is queried in
-order to determine which pixels are contained within each of the 30
-slices. Pixels within each slice are set to False (do not use) in the
-mask.
+.. note::
 
-MOS/FS Slits
-------------
-The footprints of each open MOS slitlet or fixed slit are flagged in
-a similar way as IFU slices. For MOS and FS images, the WCS object is
-queried to determine which pixels are contained within each open
-slit/slitlet and they are set to False in the mask.
+   This step, which is intended to be called in the
+   :ref:`calwebb_spec2 <calwebb_spec2>` pipeline on NIRSpec rate data,
+   is now deprecated.  In future builds, it will be replaced by
+   the :ref:`clean_flicker_noise <clean_flicker_noise_step>`
+   step, called in the :ref:`calwebb_detector1 <calwebb_detector1>`
+   pipeline on ramp data.
 
-MSA Failed Open Shutters
-------------------------
-Pixels affected by stuck open MSA shutters are masked, because they
-may contain signal. This is accomplished by setting all pixels flagged by the
-:ref:`msaflagopen <msaflagopen_step>` step with DQ value "MSA_FAILED_OPEN"
-to False in the mask.
-
-NaN Pixels
-----------
-Any pixel in the input image that has a value of NaN is temporarily reset
-to zero for input to the fitting routine and flagged as False in the mask.
-Upon completion of the noise subtraction, this population of pixels is
-set back to NaN again in the output (corrected) image.
-
-Fixed-Slit Region Pixels
-------------------------
-Full-frame MOS and IFU images may contain signal from the always open
-fixed slits, which appear in fixed region in the middle of each image.
-The entire region containing the fixed slits is masked out when
-processing MOS and IFU images. The masked region is currently hardwired
-in the step to image indexes [1:2048, 923:1116], where the indexes are
-in x, y order and in 1-indexed values.
-
-Note, however, that it is possible to plan one or more fixed slit targets
-alongside MSA slitlets in MOS observations. In this situation, the fixed
-slit region is not automatically masked.
-
-Left/Right Reference Pixel Columns
-----------------------------------
-Full-frame images contain 4 columns of reference pixels on the left and
-right edges of the image. These are not to be used in the fitting
-algorithm and hence are set to False in the mask.
-
-Outliers
---------
-Pixels in the unilluminated areas of the region can contain anomalous
-signal, due to uncaught Cosmic Rays, hot pixels, etc. A sigma-clipping
-routine is employed to find such outliers within the input image and set
-them to False in the mask. All pixels with values greater than
-:math:`median+n_sigma*sigma` are set to False in the mask.
-Here `median` and `sigma` are computed
-from the image using the astropy.stats `sigma_clipped_stats` routine,
-using the image mask to exclude pixels that have already been flagged
-and a clipping level of 5 sigma. `n_sigma` is a user-settable step
-parameter, with a default value of 5.0
-(see :ref:`nsclean step arguments <nsclean_arguments>`).
-
-Mode-Specific Masking Steps
----------------------------
-The following table indicates which flavors of masking are applied to
-images from each type of observing mode.
-
-.. |c| unicode:: U+2713 .. checkmark
-
-+--------------------------+-----+-----+-----+
-|                          |     | Mode|     |
-+--------------------------+-----+-----+-----+
-| Masking                  | IFU | MOS |  FS |
-+==========================+=====+=====+=====+
-| IFU Slices\ :sup:`1`     | |c| |     |     |
-+--------------------------+-----+-----+-----+
-| Slits/Slitlets\ :sup:`1` |     | |c| | |c| |
-+--------------------------+-----+-----+-----+
-| MSA_FAILED_OPEN          | |c| | |c| | |c| |
-+--------------------------+-----+-----+-----+
-| NaN Pixels               | |c| | |c| | |c| |
-+--------------------------+-----+-----+-----+
-| FS Region                | |c| | |c| |     |
-+--------------------------+-----+-----+-----+
-| Reference Pix            | |c| | |c| | |c| |
-+--------------------------+-----+-----+-----+
-| Outliers                 | |c| | |c| | |c| |
-+--------------------------+-----+-----+-----+
-
-:sup:`1`\ The application of these steps can be turned on and off via
-the step parameter `mask_spectral_regions`. This parameter controls
-whether the "IFU Slices" and "Slits/Slitlets" portions of the masking
-are applied.
+Overview
+========
+This step currently runs a version of the
+:ref:`clean_flicker_noise <clean_flicker_noise_step>` algorithm,
+with slightly different parameters and default values, intended
+to be backwards-compatible with the previous implementation of
+the step developed specifically for NIRSpec data. The only
+significant difference is that in this step, the cleaning is
+performed directly on the rate data, rather than group images.
+
+See the :ref:`clean_flicker_noise <clean_flicker_noise_step>`
+documentation for more information, or the
+:ref:`nsclean step arguments <nsclean_arguments>` for the default
+values used for this step.
diff --git a/docs/jwst/package_index.rst b/docs/jwst/package_index.rst
index da12f11e5..4bd3b28ba 100644
--- a/docs/jwst/package_index.rst
+++ b/docs/jwst/package_index.rst
@@ -16,6 +16,7 @@ Package Index
    badpix_selfcal/index.rst
    barshadow/index.rst
    charge_migration/index.rst
+   clean_flicker_noise/index.rst
    combine_1d/index.rst
    cube_build/index.rst
    dark_current/index.rst
diff --git a/docs/jwst/pipeline/calwebb_detector1.rst b/docs/jwst/pipeline/calwebb_detector1.rst
index e66d135af..5083b0a61 100644
--- a/docs/jwst/pipeline/calwebb_detector1.rst
+++ b/docs/jwst/pipeline/calwebb_detector1.rst
@@ -45,52 +45,52 @@ on either the first group or frame zero pixel values.
 
 .. |check| unicode:: U+2713 .. checkmark
 
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| Near-IR                                                                                  | MIRI                                                        |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| Step                                                                 | Non-TSO | TSO     | Step                                    | Non-TSO | TSO     |
-+======================================================================+=========+=========+=========================================+=========+=========+
-| :ref:`group_scale <group_scale_step>`                                | |check| | |check| | :ref:`group_scale <group_scale_step>`   | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`dq_init <dq_init_step>`                                        | |check| | |check| | :ref:`dq_init <dq_init_step>`           | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-|                                                                      |         |         | :ref:`emicorr <emicorr_step>`           | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`saturation <saturation_step>`                                  | |check| | |check| | :ref:`saturation <saturation_step>`     | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`ipc <ipc_step>` [1]_                                           |         |         | :ref:`ipc <ipc_step>`                   |         |         |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`superbias <superbias_step>`                                    | |check| | |check| | :ref:`firstframe <firstframe_step>`     | |check| |         |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`refpix <refpix_step>`                                          | |check| | |check| | :ref:`lastframe <lastframe_step>`       | |check| |         |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-|                                                                      |         |         | :ref:`reset <reset_step>`               | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`linearity <linearity_step>`                                    | |check| | |check| | :ref:`linearity <linearity_step>`       | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`persistence <persistence_step>` [2]_                           | |check| |         | :ref:`rscd <rscd_step>`                 | |check| |         |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`dark_current <dark_current_step>`                              | |check| | |check| | :ref:`dark_current <dark_current_step>` | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-|                                                                      |         |         | :ref:`refpix <refpix_step>`             | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`charge_migration <charge_migration_step>` [3]_                 | |check| |         |                                         |         |         |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`jump <jump_step>`                                              | |check| | |check| | :ref:`jump <jump_step>`                 | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`ramp_fitting <ramp_fitting_step>`                              | |check| | |check| | :ref:`ramp_fitting <ramp_fitting_step>` | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`gain_scale <gain_scale_step>`                                  | |check| | |check| | :ref:`gain_scale <gain_scale_step>`     | |check| | |check| |
-+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-
-
-.. [1] By default, the parameter reference `pars-detector1pipeline`
-   retrieved from CRDS will skip the :ref:`ipc <ipc_step>` step for all instruments.
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| Near-IR                                                                                  | MIRI                                                                     |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| Step                                                                 | Non-TSO | TSO     | Step                                                 | Non-TSO | TSO     |
++======================================================================+=========+=========+======================================================+=========+=========+
+| :ref:`group_scale <group_scale_step>`                                | |check| | |check| | :ref:`group_scale <group_scale_step>`                | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`dq_init <dq_init_step>`                                        | |check| | |check| | :ref:`dq_init <dq_init_step>`                        | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+|                                                                      |         |         | :ref:`emicorr <emicorr_step>`                        | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`saturation <saturation_step>`                                  | |check| | |check| | :ref:`saturation <saturation_step>`                  | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`ipc <ipc_step>` [1]_                                           |         |         | :ref:`ipc <ipc_step>`                                |         |         |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`superbias <superbias_step>`                                    | |check| | |check| | :ref:`firstframe <firstframe_step>`                  | |check| |         |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`refpix <refpix_step>`                                          | |check| | |check| | :ref:`lastframe <lastframe_step>`                    | |check| |         |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+|                                                                      |         |         | :ref:`reset <reset_step>`                            | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`linearity <linearity_step>`                                    | |check| | |check| | :ref:`linearity <linearity_step>`                    | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`persistence <persistence_step>` [2]_                           | |check| |         | :ref:`rscd <rscd_step>`                              | |check| |         |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`dark_current <dark_current_step>`                              | |check| | |check| | :ref:`dark_current <dark_current_step>`              | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+|                                                                      |         |         | :ref:`refpix <refpix_step>`                          | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`charge_migration <charge_migration_step>` [1]_                 | |check| |         |                                                      |         |         |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`jump <jump_step>`                                              | |check| | |check| | :ref:`jump <jump_step>`                              | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`clean_flicker_noise <clean_flicker_noise_step>` [1]_           | |check| | |check| | :ref:`clean_flicker_noise <clean_flicker_noise_step>`| |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`ramp_fitting <ramp_fitting_step>`                              | |check| | |check| | :ref:`ramp_fitting <ramp_fitting_step>`              | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+| :ref:`gain_scale <gain_scale_step>`                                  | |check| | |check| | :ref:`gain_scale <gain_scale_step>`                  | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+
+
+
+.. [1] By default, this step is skipped in the ``calwebb_detector1`` pipeline
+   for all instruments and modes.
 .. [2] The :ref:`persistence <persistence_step>` step is currently hardwired to be skipped in
    the `Detector1Pipeline` module for all NIRSpec exposures.
-.. [3] By default, the :ref:`charge_migration <charge_migration_step>` step is skipped in
-   the `Detector1Pipeline` module for all instruments.
-   
+
 
 Arguments
 ---------
diff --git a/jwst/clean_flicker_noise/__init__.py b/jwst/clean_flicker_noise/__init__.py
new file mode 100644
index 000000000..e9bd9a65d
--- /dev/null
+++ b/jwst/clean_flicker_noise/__init__.py
@@ -0,0 +1,3 @@
+from .clean_flicker_noise_step import CleanFlickerNoiseStep
+
+__all__ = ['CleanFlickerNoiseStep']
diff --git a/jwst/clean_flicker_noise/clean_flicker_noise.py b/jwst/clean_flicker_noise/clean_flicker_noise.py
new file mode 100644
index 000000000..2a41da373
--- /dev/null
+++ b/jwst/clean_flicker_noise/clean_flicker_noise.py
@@ -0,0 +1,1522 @@
+import logging
+import warnings
+
+import gwcs
+import numpy as np
+from astropy.stats import sigma_clipped_stats, SigmaClip
+from astropy.utils.exceptions import AstropyUserWarning
+from photutils.background import Background2D, MedianBackground
+from scipy.optimize import curve_fit
+from stdatamodels.jwst.datamodels import dqflags
+
+from jwst import datamodels
+from jwst.assign_wcs import nirspec, AssignWcsStep
+from jwst.flatfield import FlatFieldStep
+from jwst.clean_flicker_noise.lib import NSClean, NSCleanSubarray
+from jwst.lib.basic_utils import LoggingContext
+from jwst.msaflagopen import MSAFlagOpenStep
+from jwst.ramp_fitting import RampFitStep
+
+
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
+# Fixed slit region to mask, for NIRSpec MOS and IFU data
+# Values are y start and stop indices, for the edges of the
+# region to mask.
+NRS_FS_REGION = [922, 1116]
+
+
+def make_rate(input_model, input_dir='', return_cube=False):
+    """
+    Make a rate model from a ramp model.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodel.RampModel`
+        Input ramp model.
+
+    input_dir : str
+        Path to the input directory.  Used by sub-steps (e.g. assign_wcs
+        for NIRSpec MOS data) to find auxiliary data.
+
+    return_cube : bool, optional
+        If set, a CubeModel will be returned, with a separate
+        rate for each integration.  Otherwise, an ImageModel is returned
+        with the combined rate for the full observation.
+
+    Returns
+    -------
+    rate_model : `~jwst.datamodel.ImageModel` or `~jwst.datamodel.CubeModel`
+        The rate or rateints model.
+    """
+    # Call the ramp fit step on a copy of the input
+    # Use software default values for parameters
+
+    log.info("Creating draft rate file for scene masking")
+    step = RampFitStep()
+    step.input_dir = input_dir
+    with LoggingContext(step.log, level=logging.WARNING):
+        # Note: the copy is currently needed because ramp fit
+        # closes the input model when it's done, and we need
+        # it to stay open.
+        rate, rateints = step.run(input_model.copy())
+
+    if return_cube:
+        output_model = rateints
+        rate.close()
+        del rate
+    else:
+        output_model = rate
+        rateints.close()
+        del rateints
+
+    return output_model
+
+
+def post_process_rate(input_model, input_dir='', assign_wcs=False,
+                      msaflagopen=False, flat_dq=False):
+    """
+    Post-process the input rate model, as needed.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodel.ImageModel` or `~jwst.datamodel.CubeModel`
+        Input rate model.
+
+    input_dir : str
+        Path to the input directory.  Used by sub-steps (e.g. assign_wcs
+        for NIRSpec MOS data) to find auxiliary data.
+
+    assign_wcs : bool, optional
+        If set and the input does not already have a WCS assigned,
+        the assign_wcs step will be called on the rate model.
+
+    msaflagopen : bool, optional
+        If set, the msaflagopen step will be called on the rate model.
+        If a WCS is not already present, assign_wcs will be called first.
+
+    flat_dq : bool, optional
+        If set, the flat_field step will be run on the input model. DQ
+        flags are retrieved from the output and added to the input
+        model's DQ array. The rate data is not modified.
+
+    Returns
+    -------
+    output_model : `~jwst.datamodel.ImageModel` or `~jwst.datamodel.CubeModel`
+        The updated model.
+    """
+    output_model = input_model
+
+    # If needed, assign a WCS
+    if (assign_wcs or msaflagopen) and not hasattr(output_model.meta, 'wcs'):
+        log.info("Assigning a WCS for scene masking")
+        step = AssignWcsStep()
+        step.input_dir = input_dir
+        with LoggingContext(step.log, level=logging.WARNING):
+            output_model = step.run(output_model)
+
+    # If needed, flag open MSA shutters
+    if msaflagopen:
+        log.info("Flagging failed-open MSA shutters for scene masking")
+        step = MSAFlagOpenStep()
+        step.input_dir = input_dir
+        with LoggingContext(step.log, level=logging.WARNING):
+            output_model = step.run(output_model)
+
+    # If needed, draft a flat correction to retrieve non-science areas
+    if flat_dq:
+        log.info("Retrieving flat DQ values for scene masking")
+        step = FlatFieldStep()
+        step.input_dir = input_dir
+        with LoggingContext(step.log, level=logging.WARNING):
+            flat_corrected_model = step.run(output_model)
+
+        # Copy out the flat DQ plane, leave the data as is
+        output_model.dq |= flat_corrected_model.dq
+        flat_corrected_model.close()
+        del flat_corrected_model
+
+    return output_model
+
+
+def mask_ifu_slices(input_model, mask):
+    """
+    Flag pixels within IFU slices.
+
+    Find pixels located within IFU slices, according to the WCS,
+    and flag them in the mask, so that they do not get used.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodel.JwstDataModel`
+        Science data model
+
+    mask : array-like of bool
+        2D input mask that will be updated. True indicates background
+        pixels to be used. IFU slice regions will be set to False.
+
+    Returns
+    -------
+    mask : array-like of bool
+        2D output mask with additional flags for science pixels
+    """
+    log.info("Finding slice pixels for an IFU image")
+
+    # Initialize global DQ map to all zero (OK to use)
+    dqmap = np.zeros_like(input_model.dq)
+
+    # Get the wcs objects for all IFU slices
+    list_of_wcs = nirspec.nrs_ifu_wcs(input_model)
+
+    # Loop over the IFU slices, finding the valid region for each
+    for (k, ifu_wcs) in enumerate(list_of_wcs):
+
+        # Construct array indexes for pixels in this slice
+        x, y = gwcs.wcstools.grid_from_bounding_box(
+            ifu_wcs.bounding_box, step=(1, 1), center=True)
+
+        # Get the world coords for all pixels in this slice;
+        # all we actually need are wavelengths
+        coords = ifu_wcs(x, y)
+        dq = dqmap[..., y.astype(int), x.astype(int)]
+        wl = coords[2]
+
+        # set non-NaN wavelength locations as do not use (one)
+        valid = ~np.isnan(wl)
+        dq = dq[..., valid]
+        x = x[..., valid]
+        y = y[..., valid]
+        dq[:] = 1
+
+        # Copy DQ for this slice into global DQ map
+        dqmap[..., y.astype(int), x.astype(int)] = dq
+
+    # Now set all non-zero locations in the mask to False (do not use)
+    mask[dqmap == 1] = False
+
+    return mask
+
+
+def mask_slits(input_model, mask):
+    """
+    Flag pixels within science regions.
+
+    Find pixels located within MOS or fixed slit footprints
+    and flag them in the mask, so that they do not get used.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodel.JwstDataModel`
+        Science data model.
+
+    mask : array-like of bool
+        2D input mask that will be updated. True indicates background
+        pixels to be used. Slit regions will be set to False.
+
+    Returns
+    -------
+    mask : array-like of bool
+        2D output mask with additional flags for slit pixels
+    """
+
+    from jwst.extract_2d.nirspec import offset_wcs
+
+    log.info("Finding slit/slitlet pixels")
+
+    # Get the slit-to-msa frame transform from the WCS object
+    slit2msa = input_model.meta.wcs.get_transform('slit_frame', 'msa_frame')
+
+    # Loop over the slits, marking all the pixels within each bounding
+    # box as False (do not use) in the mask.
+    # Note that for 3D masks (TSO mode), all planes will be set to the same value.
+    for slit in slit2msa.slits:
+        slit_wcs = nirspec.nrs_wcs_set_input(input_model, slit.name)
+        xlo, xhi, ylo, yhi = offset_wcs(slit_wcs)
+        mask[..., ylo:yhi, xlo:xhi] = False
+
+    return mask
+
+
+def clip_to_background(image, mask, sigma_lower=3.0, sigma_upper=2.0,
+                       fit_histogram=False, lower_half_only=False,
+                       verbose=False):
+    """
+    Flag signal and bad pixels in the image mask.
+
+    Given an image, estimate the background level and a sigma value for the
+    mean background.
+
+    The center and sigma may be calculated from a simple sigma-clipped
+    median and standard deviation, or may be refined by fitting a Gaussian
+    to a histogram of the values.  In that case, the center is the
+    Gaussian mean and the sigma is the Gaussian width.
+
+    Pixels above the center + sigma_upper * sigma are assumed to
+    have signal; those below this level are assumed to be background pixels.
+
+    Pixels less than center - sigma_lower * sigma are also excluded as bad values.
+
+    The input mask is updated in place.
+
+    Parameters
+    ----------
+    image : array-like of float
+        2D image containing signal and background values.
+
+    mask : array-like of bool
+        2D input mask to be updated. True indicates background
+        pixels to be used. Regions containing signal or outliers
+        will be set to False.
+
+    sigma_lower : float, optional
+        The number of standard deviations to use as the lower bound
+        for the clipping limit. Values below this limit are marked
+        False in the mask.
+
+    sigma_upper : float, optional
+        The number of standard deviations to use as the upper bound
+        for the clipping limit. Values above this limit are marked
+        False in the mask.
+
+    fit_histogram :  bool, optional
+        If set, the center value and standard deviation used with
+        `sigma_lower` and `sigma_upper` for clipping outliers is derived
+        from a Gaussian fit to a histogram of values. Otherwise, the
+        center and standard deviation are derived from a simple iterative
+        sigma clipping.
+
+    lower_half_only : bool, optional
+        If set, the data used to compute the center and standard deviation
+        for clipping is the lower half of the distribution only. Values
+        below the median are mirrored around the median value to simulate
+        a symmetric distribution.  This is intended to account for
+        asymmetrical value distributions, with long tails in the upper
+        half of the distribution, due to diffuse emission, for example.
+
+    verbose : bool, optional
+        If set, DEBUG level messages are issued with details on the
+        computed statistics.
+    """
+    # Sigma limit for basic stats
+    sigma_limit = 3.0
+
+    # Check mask for any valid data before proceeding
+    if not np.any(mask):
+        return
+
+    # Initial iterative sigma clip
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            action="ignore", category=AstropyUserWarning)
+        mean, median, sigma = sigma_clipped_stats(
+            image, mask=~mask, sigma=sigma_limit)
+    if fit_histogram:
+        center = mean
+    else:
+        center = median
+    if verbose:
+        log.debug('From initial sigma clip:')
+        log.debug(f'    center: {center:.5g}')
+        log.debug(f'    sigma: {sigma:.5g}')
+
+    # If desired, use only the lower half of the data distribution
+    if lower_half_only:
+        lower_half_idx = mask & (image <= center)
+        data_for_stats = np.concatenate(
+            ((image[lower_half_idx] - center),
+             (center - image[lower_half_idx]))) + center
+
+        # Redo stats on lower half of distribution
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                action="ignore", category=AstropyUserWarning)
+            mean, median, sigma = sigma_clipped_stats(
+                data_for_stats, sigma=sigma_limit)
+        if fit_histogram:
+            center = mean
+        else:
+            center = median
+        if verbose:
+            log.debug('From lower half distribution:')
+            log.debug(f'    center: {center:.5g}')
+            log.debug(f'    sigma: {sigma:.5g}')
+    else:
+        data_for_stats = image[mask]
+
+    # Refine sigma and center from a fit to a histogram, if desired
+    if fit_histogram:
+        try:
+            hist, edges = np.histogram(
+                data_for_stats, bins=2000,
+                range=(center - 4. * sigma, center + 4. * sigma))
+        except ValueError:
+            log.error('Histogram failed; using clip center and sigma.')
+            hist, edges = None, None
+
+        param_opt = None
+        if hist is not None:
+            values = (edges[1:] + edges[0:-1]) / 2.
+            ind = np.argmax(hist)
+            mode_estimate = values[ind]
+
+            # Fit a Gaussian profile to the histogram
+            def gaussian(x, g_amp, g_mean, g_sigma):
+                return g_amp * np.exp(-0.5 * ((x - g_mean) / g_sigma) ** 2)
+
+            param_start = (hist[ind], mode_estimate, sigma)
+            bounds = [(0, values[0], 0),
+                      (np.inf, values[-1], values[-1] - values[0])]
+            try:
+                param_opt, _ = curve_fit(gaussian, values, hist, p0=param_start,
+                                         bounds=bounds)
+            except RuntimeError:
+                log.error('Gaussian fit failed; using clip center and sigma.')
+                param_opt = None
+
+            if verbose:
+                log.debug('From histogram:')
+                log.debug(f'    mode estimate: {mode_estimate:.5g}')
+                log.debug(f'    range of values in histogram: '
+                          f'{values[0]:.5g} to {values[-1]:.5g}')
+
+        if verbose:
+            log.debug('Gaussian fit results:')
+        if param_opt is None:
+            if verbose:
+                log.debug('    (fit failed)')
+        else:
+            if verbose:
+                log.debug(f'    peak: {param_opt[0]:.5g}')
+                log.debug(f'    center: {param_opt[1]:.5g}')
+                log.debug(f'    sigma: {param_opt[2]:.5g}')
+            center = param_opt[1]
+            sigma = param_opt[2]
+
+    # Set limits from center and sigma
+    background_lower_limit = center - sigma_lower * sigma
+    background_upper_limit = center + sigma_upper * sigma
+    if verbose:
+        log.debug(f'Mask limits: {background_lower_limit:.5g} '
+                  f'to {background_upper_limit:.5g}')
+
+    # Clip bad values
+    bad_values = image < background_lower_limit
+    mask[bad_values] = False
+
+    # Clip signal (> N sigma)
+    signal = image > background_upper_limit
+    mask[signal] = False
+
+
+def create_mask(input_model, mask_science_regions=False,
+                n_sigma=2.0, fit_histogram=False, single_mask=False):
+    """
+    Create a mask identifying background pixels.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodel.JwstDataModel`
+        Science data model, containing rate data with all necessary
+        pre-processing already performed.
+
+    mask_science_regions : bool, optional
+        For NIRSpec, mask regions of the image defined by WCS bounding
+        boxes for slits/slices, as well as any regions known to be
+        affected by failed-open MSA shutters. This requires the
+        `assign_wcs` and `msaflagopen` steps to have been run on the
+        input_model.
+        For MIRI imaging, mask regions of the detector not used for science.
+        This requires that DO_NOT_USE flags are set in the DQ array
+        for the input_model.
+
+    n_sigma : float, optional
+        Sigma threshold for masking outliers.
+
+    fit_histogram : bool, optional
+        If set, the 'sigma' used with `n_sigma` for clipping outliers
+        is derived from a Gaussian fit to a histogram of values.
+        Otherwise, a simple iterative sigma clipping is performed.
+
+    single_mask : bool, optional
+        If set, a single mask will be created, regardless of
+        the number of input integrations. Otherwise, the mask will
+        be a 3D cube, with one plane for each integration.
+
+    Returns
+    -------
+    mask : array-like of bool
+        2D or 3D image mask
+    """
+    exptype = input_model.meta.exposure.type.lower()
+
+    # Initialize mask to all True. Subsequent operations will mask
+    # out pixels that contain signal.
+    mask = np.full(input_model.dq.shape, True)
+
+    # If IFU, mask all pixels contained in the IFU slices
+    if exptype == 'nrs_ifu' and mask_science_regions:
+        mask = mask_ifu_slices(input_model, mask)
+
+    # If MOS or FS, mask all pixels affected by open slitlets
+    if (exptype in ['nrs_fixedslit', 'nrs_brightobj', 'nrs_msaspec']
+            and mask_science_regions):
+        mask = mask_slits(input_model, mask)
+
+    # If IFU or MOS, mask pixels affected by failed-open shutters
+    if mask_science_regions and exptype in ['nrs_ifu', 'nrs_msaspec']:
+        open_pix = input_model.dq & dqflags.pixel['MSA_FAILED_OPEN']
+        mask[open_pix > 0] = False
+
+    # If MIRI imaging, mask the non-science regions:
+    # they contain irrelevant emission
+    if mask_science_regions and exptype in ['mir_image']:
+        non_science = (input_model.dq & dqflags.pixel['DO_NOT_USE']) > 0
+        mask[non_science] = False
+
+    # Mask any NaN pixels or exactly zero value pixels
+    no_data_pix = np.isnan(input_model.data) | (input_model.data == 0)
+    mask[no_data_pix] = False
+
+    # If IFU or MOS, mask the fixed-slit area of the image; uses hardwired indexes
+    if mask_science_regions:
+        if exptype == 'nrs_ifu':
+            log.info("Masking the fixed slit region for IFU data.")
+            mask[..., NRS_FS_REGION[0]:NRS_FS_REGION[1], :] = False
+        elif exptype == 'nrs_msaspec':
+            # check for any slits defined in the fixed slit quadrant:
+            # if there is nothing there of interest, mask the whole FS region
+            try:
+                slit2msa = input_model.meta.wcs.get_transform('slit_frame', 'msa_frame')
+                is_fs = [s.quadrant == 5 for s in slit2msa.slits]
+            except (AttributeError, ValueError, TypeError):
+                log.warning('Slit to MSA transform not found.')
+                is_fs = [False]
+            if not any(is_fs):
+                log.info("Masking the fixed slit region for MOS data.")
+                mask[..., NRS_FS_REGION[0]:NRS_FS_REGION[1], :] = False
+            else:
+                log.info("Fixed slits found in MSA definition; "
+                         "not masking the fixed slit region for MOS data.")
+
+    # Mask outliers and signal using sigma clipping stats.
+    # For 3D data, loop over each integration separately.
+    if input_model.data.ndim == 3:
+        for i in range(input_model.data.shape[0]):
+            clip_to_background(
+                input_model.data[i], mask[i],
+                sigma_upper=n_sigma, fit_histogram=fit_histogram, verbose=True)
+    else:
+        clip_to_background(
+            input_model.data, mask,
+            sigma_upper=n_sigma, fit_histogram=fit_histogram, verbose=True)
+
+    # Reduce the mask to a single plane if needed
+    if single_mask and mask.ndim == 3:
+        mask = np.all(mask, axis=0)
+
+    # Return the mask
+    return mask
+
+
+def background_level(image, mask, background_method='median',
+                     background_box_size=None):
+    """
+    Fit a low-resolution background level.
+
+    Parameters
+    ----------
+    image : array-like of float
+        The 2D image containing the background to fit.
+
+    mask : array-like of bool
+        The mask that indicates which pixels are to be used in fitting.
+        True indicates a background pixel.
+
+    background_method : {'median', 'model', None}, optional
+        If 'median', the preliminary background to remove and restore
+        is a simple median of the background data.  If 'model', the
+        background data is fit with a low-resolution model via
+        `~photutils.background.Background2D`.  If None, the background
+        value is 0.0.
+
+    background_box_size : tuple of int, optional
+        Box size for the data grid used by `Background2D` when
+        `background_method` = 'model'. For best results, use a box size
+        that evenly divides the input image shape. Defaults to 32x32 if
+        not provided.
+
+    Returns
+    -------
+    background : float or array-like of float
+        The background level: a single value, if `background_method`
+        is 'median' or None, or an array matching the input image size
+        if `background_method` is 'model'.
+    """
+    if background_method is None:
+        background = 0.0
+
+    else:
+        # Sigma limit for basic stats
+        sigma_limit = 3.0
+
+        # Flag more signal in the background subtracted image,
+        # with sigma set by the lower half of the distribution only
+        clip_to_background(
+            image, mask, sigma_lower=sigma_limit, sigma_upper=sigma_limit,
+            lower_half_only=True)
+
+        if background_method == 'model':
+            sigma_clip_for_bkg = SigmaClip(sigma=sigma_limit, maxiters=5)
+            bkg_estimator = MedianBackground()
+
+            if background_box_size is None:
+                background_box_size = [32, 32]
+            box_division_remainder = (image.shape[0] % background_box_size[0],
+                                      image.shape[1] % background_box_size[1])
+            if not np.allclose(box_division_remainder, 0):
+                log.warning(f'Background box size {background_box_size} '
+                            f'does not divide evenly into the image '
+                            f'shape {image.shape}.')
+
+            try:
+                with warnings.catch_warnings():
+                    warnings.filterwarnings(
+                        action="ignore", category=AstropyUserWarning)
+                    bkg = Background2D(
+                        image, box_size=background_box_size, filter_size=(5, 5),
+                        mask=~mask, sigma_clip=sigma_clip_for_bkg,
+                        bkg_estimator=bkg_estimator)
+                background = bkg.background
+            except ValueError:
+                log.error('Background fit failed, using median value.')
+                background = np.nanmedian(image[mask])
+        else:
+            background = np.nanmedian(image[mask])
+    return background
+
+
+def fft_clean_full_frame(image, mask, detector):
+    """
+    Fit and remove background noise in frequency space for a full-frame image.
+
+    Parameters
+    ----------
+    image : array-like of float
+        The image to be cleaned.
+
+    mask : array-like of bool
+        The mask that indicates which pixels are to be used in fitting.
+
+    detector : str
+        The name of the detector from which the data originate.
+
+    Returns
+    -------
+    cleaned_image : array-like of float
+        The cleaned image.
+    """
+
+    # Instantiate the cleaner
+    cleaner = NSClean(detector, mask)
+
+    # Clean the image
+    try:
+        cleaned_image = cleaner.clean(image, buff=True)
+    except np.linalg.LinAlgError:
+        log.warning("Error cleaning image; step will be skipped")
+        return None
+
+    return cleaned_image
+
+
+def fft_clean_subarray(image, mask, detector, npix_iter=512,
+                       fc=(1061, 1211, 49943, 49957),
+                       exclude_outliers=True, sigrej=4, minfrac=0.05):
+    """
+    Fit and remove background noise in frequency space for a subarray image.
+
+    Parameters
+    ----------
+    image : array-like of float
+        The image to be cleaned.
+
+    mask : array-like of bool
+        The mask that indicates which pixels are to be used in fitting.
+
+    detector : str
+        The name of the detector from which the data originate.
+
+    npix_iter : int
+        Number of pixels to process simultaneously.  Default 512.  Should
+        be at least a few hundred to access sub-kHz frequencies in areas
+        where most pixels are available for fitting.  Previous default
+        behavior corresponds to npix_iter of infinity.
+
+    fc : tuple
+        Apodizing filter definition. These parameters are tunable. The
+        defaults happen to work well for NIRSpec BOTS exposures.
+          1) Unity gain for f < fc[0]
+          2) Cosine roll-off from fc[0] to fc[1]
+          3) Zero gain from fc[1] to fc[2]
+          4) Cosine roll-on from fc[2] to fc[3]
+        Default (1061, 1211, 49943, 49957)
+
+    exclude_outliers : bool
+        Find and mask outliers in the fit?  Default True
+
+    sigrej : float
+        Number of sigma to clip when identifying outliers.  Default 4.
+
+    minfrac : float
+        Minimum fraction of pixels locally available in the mask in
+        order to attempt a correction.  Default 0.05 (i.e., 5%).
+
+    Returns
+    -------
+    cleaned_image : array-like of float
+        The cleaned image.
+    """
+    # Flip the image to detector coords. NRS1 requires a transpose
+    # of the axes, while NRS2 requires a transpose and flip.
+    if detector == "NRS1":
+        image = image.transpose()
+        mask = mask.transpose()
+    else:
+        image = image.transpose()[::-1]
+        mask = mask.transpose()[::-1]
+
+    # We must do the masking of discrepant pixels here: it just
+    # doesn't work if we wait and do it in the cleaner.  This is
+    # basically copied from lib.py.  Use a robust estimator for
+    # standard deviation, then exclude discrepant pixels and their
+    # four nearest neighbors from the fit.
+    
+    if exclude_outliers:
+        med = np.median(image[mask])
+        std = 1.4825 * np.median(np.abs((image - med)[mask]))
+        outlier = mask & (np.abs(image - med) > sigrej * std)
+        
+        mask = mask & (~outlier)
+        
+        # also get four nearest neighbors of flagged pixels
+        mask[1:] = mask[1:] & (~outlier[:-1])
+        mask[:-1] = mask[:-1] & (~outlier[1:])
+        mask[:, 1:] = mask[:, 1:] & (~outlier[:, :-1])
+        mask[:, :-1] = mask[:, :-1] & (~outlier[:, 1:])
+
+    # Used to determine the fitting intervals along the slow scan
+    # direction.  Pre-pend a zero so that sum_mask[i] is equal
+    # to np.sum(mask[:i], axis=1).
+
+    sum_mask = np.array([0] + list(np.cumsum(np.sum(mask, axis=1))))
+
+    # i1 will be the first row with a nonzero element in the mask
+    # imax will be the last row with a nonzero element in the mask
+    
+    nonzero_mask_element = np.sum(mask, axis=1) > 0
+
+    if np.sum(nonzero_mask_element) == 0:
+        log.warning("No good pixels in mask; step will be skipped")
+        return None
+
+    i1 = np.amin(np.arange(mask.shape[0])[nonzero_mask_element])
+    imax = np.amax(np.arange(mask.shape[0])[nonzero_mask_element])
+        
+    i1_vals = []
+    di_list = []
+    models = []
+    while i1 <= imax:
+
+        # Want npix_iter available pixels in this section.  If
+        # there are fewer than 1.5*npix_iter available pixels in
+        # the rest of the image, just go to the end.
+        k = 0
+        for k in range(i1 + 1, imax + 2):
+            if (sum_mask[k] - sum_mask[i1] > npix_iter
+                    and sum_mask[-1] - sum_mask[i1] > 1.5 * npix_iter):
+                break
+            
+        di = k - i1
+
+        i1_vals += [i1]
+        di_list += [di]
+
+        # Fit this section only if at least minpct% of the pixels
+        # are available for finding the background.  Don't flag
+        # outliers section-by-section; we have to do that earlier
+        # over the full array to get reliable values for the mean
+        # and standard deviation.
+        
+        if np.mean(mask[i1:i1 + di]) > minfrac:
+            cleaner = NSCleanSubarray(image[i1:i1 + di], mask[i1:i1 + di],
+                                      fc=fc, exclude_outliers=False)
+            try:
+                models += [cleaner.clean(return_model=True)]
+            except np.linalg.LinAlgError:
+                log.warning("Error cleaning image; step will be skipped")
+                return None
+        else:
+            log.warning("Insufficient reference pixels for NSClean around "
+                        "row %d; no correction will be made here." % i1)
+            models += [np.zeros(image[i1:i1 + di].shape)]
+
+        # If we have reached the end of the array, we are finished.
+        if k == imax + 1:
+            break
+
+        # Step forward by half an interval so that we have
+        # overlapping fitting regions.
+        
+        i1 += max(int(np.round(di/2)), 1)
+            
+    model = np.zeros(image.shape)
+    tot_wgt = np.zeros(image.shape)
+
+    # When we combine different corrections computed over
+    # different intervals, each one the highest weight towards the
+    # center of its interval and less weight towards the edge.
+    # Use nonzero weights everywhere so that if only one
+    # correction is available it gets unit weight when we
+    # normalize.
+    
+    for i in range(len(models)):
+        wgt = 1.001 - np.abs(np.linspace(-1, 1, di_list[i]))[:, np.newaxis]
+        model[i1_vals[i]:i1_vals[i] + di_list[i]] += wgt*models[i]
+        tot_wgt[i1_vals[i]:i1_vals[i] + di_list[i]] += wgt
+
+    # don't divide by zero
+    tot_wgt[model == 0] = 1
+    model /= tot_wgt
+    cleaned_image = image - model
+
+    # Restore the cleaned image to the science frame
+    if detector == "NRS1":
+        cleaned_image = cleaned_image.transpose()
+    else:
+        cleaned_image = cleaned_image[::-1].transpose()
+
+    return cleaned_image
+
+
+def median_clean(image, mask, axis_to_correct, fit_by_channel=False):
+    """
+    Fit and remove background noise via median values along one image axis.
+
+    Parameters
+    ----------
+    image : array-like of float
+        The image to be cleaned.
+
+    mask : array-like of bool
+        The mask that indicates which pixels are to be used in fitting.
+        True indicates a background pixel.
+
+    axis_to_correct : int
+        For NIR detectors, the axis to correct should be the detector slow
+        readout direction.  Values expected are 1 or 2, following the JWST
+        datamodel definition (meta.subarray.slowaxis).  For MIRI, flicker
+        noise appears along the vertical direction, so `axis_to_correct`
+        should be set to 1 (median along the y-axis).
+
+    fit_by_channel : bool, optional
+        If set, flicker noise is fit independently for each detector channel.
+        Ignored for MIRI and for subarray data.
+
+    Returns
+    -------
+    cleaned_image : array-like of float
+        The cleaned image.
+    """
+    # Masked median along slow axis
+    masked_image = np.ma.array(image, mask=~mask)
+    corrected_image = image.copy()
+    array_axis = axis_to_correct - 1
+
+    # If desired, take the median over each channel separately.
+
+    # Full frame for the NIR detectors is 2048 x 2048, and there
+    # are 4 channels per detector, divided along the slow axis.
+    # The fit_by_channel option should only be set for full frame
+    # NIR data, but make sure the slow axis has the right size
+    # here anyway.
+    if fit_by_channel and image.shape[array_axis] == 2048:
+        n_output = 4
+        channel_size = 512
+    else:
+        # For MIRI data and subarrays, the whole image should
+        # be considered at once, as a single "channel"
+        n_output = 1
+        channel_size = image.shape[array_axis]
+
+    # Compute stripes from the median over the background data
+    # in the channel
+    cstart = 0
+    cstop = channel_size
+    for channel in range(n_output):
+        if array_axis == 1:
+            channel_image = masked_image[:, cstart:cstop]
+        else:
+            channel_image = masked_image[cstart:cstop, :]
+        stripes = np.ma.median(channel_image, axis=array_axis, keepdims=True)
+        stripes = np.ma.filled(stripes, fill_value=0.0)
+
+        # Remove median stripes
+        if array_axis == 1:
+            corrected_image[:, cstart:cstop] = image[:, cstart:cstop] - stripes
+        else:
+            corrected_image[cstart:cstop, :] = image[cstart:cstop, :] - stripes
+
+        cstart += channel_size
+        cstop += channel_size
+
+    return corrected_image
+
+
+def _check_input(exp_type, fit_method):
+    """
+    Check for valid input data and options.
+
+    Parameters
+    ----------
+    exp_type : str
+        Exposure type for the input.
+    fit_method : str
+        Noise fitting method.
+
+    Returns
+    -------
+    valid : bool
+        True if the input is valid.
+    """
+    message = None
+    miri_allowed = ['MIR_IMAGE']
+    if exp_type.startswith('MIR') and exp_type not in miri_allowed:
+        message = f"EXP_TYPE {exp_type} is not supported"
+
+    nsclean_allowed = ['NRS_MSASPEC', 'NRS_IFU', 'NRS_FIXEDSLIT', 'NRS_BRIGHTOBJ']
+    if fit_method == 'fft':
+        if exp_type not in nsclean_allowed:
+            message = f"Fit method 'fft' cannot be applied to exp_type {exp_type}"
+
+    if message is not None:
+        log.warning(message)
+        log.warning("Step will be skipped")
+        return False
+
+    return True
+
+
+def _make_intermediate_model(input_model, intermediate_data):
+    """
+    Make a data model to contain intermediate outputs.
+
+    The output model type depends on the shape of the input
+    intermediate data.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodel.JwstDataModel`
+        The input data.
+    intermediate_data : array-like
+        The intermediate data to save.
+
+    Returns
+    -------
+    intermediate_model : ~jwst.datamodel.JwstDataModel`
+        A model containing only the intermediate data and top-level
+        metadata matching the input.
+    """
+    if intermediate_data.ndim == 4:
+        intermediate_model = datamodels.RampModel(data=intermediate_data)
+    elif intermediate_data.ndim == 3:
+        intermediate_model = datamodels.CubeModel(data=intermediate_data)
+    else:
+        intermediate_model = datamodels.ImageModel(data=intermediate_data)
+
+    # Copy metadata from input model
+    intermediate_model.update(input_model)
+    return intermediate_model
+
+
+def _standardize_parameters(
+        exp_type, subarray, slowaxis, background_method, fit_by_channel):
+    """
+    Standardize input parameters.
+
+    Check input parameters against input exposure type and assemble
+    values needed for subsequent correction.
+
+    Parameters
+    ----------
+    exp_type : str
+        Exposure type for the input data.
+    subarray : str
+        Subarray name for the input data.
+    slowaxis : int
+        Detector slow axis.
+    background_method : str
+        Input option for background method.
+    fit_by_channel : bool
+        Input option to fit noise by channel.
+
+    Returns
+    -------
+    axis_to_correct : int
+        Axis along which flicker noise appears for the input
+        exposure type.
+    background_method : str or None
+        Standardized parameter value for background correction.
+    fit_by_channel : bool
+        Standardized parameter value for the input subarray and
+        exposure type.
+    fc : tuple of int
+        Frequency cutoff values for use with FFT correction,
+        by input subarray.
+    """
+    # Get axis to correct, by instrument
+    if exp_type.startswith('MIR'):
+        # MIRI doesn't have 1/f-noise, but it does have a vertical flickering.
+        # Set the axis for median correction to the y-axis.
+        axis_to_correct = 1
+    else:
+        # For NIR detectors, the axis for median correction is the
+        # detector slowaxis.
+        axis_to_correct = abs(slowaxis)
+
+    # Standardize background arguments
+    if str(background_method).lower() == 'none':
+        background_method = None
+
+    # Check for fit_by_channel argument, and use only if data is full frame
+    if fit_by_channel and (subarray != 'FULL' or exp_type.startswith('MIR')):
+        log.warning("Fit by channel can only be used for full-frame NIR data.")
+        log.warning("Setting fit_by_channel to False.")
+        fit_by_channel = False
+
+    # For a fft correction, ALLSLITS exposures need different
+    # ranges of 1/f frequencies.  Be less aggressive with
+    # fitting higher frequencies in ALLSLITS mode.
+    if subarray == "ALLSLITS":
+        fc = (150, 200, 49943, 49957)
+    else:
+        fc = (1061, 1211, 49943, 49957)
+
+    return axis_to_correct, background_method, fit_by_channel, fc
+
+
+def _make_processed_rate_image(
+        input_model, single_mask, input_dir, exp_type, mask_science_regions):
+    """
+    Make a draft rate image and postprocess if needed.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodel.JwstDataModel`
+        The input data.
+    single_mask : bool
+        If set, a single scene mask is desired, so create
+        a single rate image.  Otherwise, create a rateints
+        cube, with independent rate information for each
+        integration.
+    input_dir : str
+        Path to the input directory.  Used by sub-steps to find auxiliary
+        data.
+    exp_type : str
+        Exposure type for the input data, used to determine
+        which kinds of postprocessing is necessary.
+    mask_science_regions : bool
+        If set, science regions should be masked, so run `assign_wcs`
+        and `msaflagopen` for NIRSpec and `flatfield` for MIRI
+        after creating the draft rate image.
+
+    Returns
+    -------
+    image_model : `~jwst.datamodel.JwstDataModel`
+        The processed rate image or cube.
+    """
+    if isinstance(input_model, datamodels.RampModel):
+        image_model = make_rate(input_model, return_cube=(not single_mask),
+                                input_dir=input_dir)
+    else:
+        # input is already a rate file
+        image_model = input_model
+
+    # If needed, assign a WCS to the rate file,
+    # flag open MSA shutters, or retrieve flat DQ flags
+    assign_wcs = exp_type.startswith('NRS')
+    flag_open = (exp_type in ['NRS_IFU', 'NRS_MSASPEC'])
+    flat_dq = (exp_type in ['MIR_IMAGE'])
+    if mask_science_regions:
+        image_model = post_process_rate(
+            image_model, assign_wcs=assign_wcs,
+            msaflagopen=flag_open, flat_dq=flat_dq,
+            input_dir=input_dir)
+
+    return image_model
+
+
+def _make_scene_mask(
+        user_mask, image_model, mask_science_regions,
+        n_sigma, fit_histogram, single_mask, save_mask):
+    """
+    Make a scene mask from user input or rate image.
+
+    If provided, the user mask is opened as a datamodel and directly
+    returned. Otherwise, the mask is generated from the rate data in
+    `image_model`.
+
+    Parameters
+    ----------
+    user_mask : str or None
+        Path to user-supplied mask image.
+    image_model : `~jwst.datamodel.JwstDataModel`
+        A rate image or cube, processed as needed.
+    mask_science_regions: bool
+        For NIRSpec, mask regions of the image defined by WCS bounding
+        boxes for slits/slices, as well as any regions known to be
+        affected by failed-open MSA shutters.  For MIRI imaging, mask
+        regions of the detector not used for science.
+    n_sigma : float
+        N-sigma rejection level for finding outliers.
+    fit_histogram : bool
+        If set, the 'sigma' used with `n_sigma` for clipping outliers
+        is derived from a Gaussian fit to a histogram of values.
+        Otherwise, a simple iterative sigma clipping is performed.
+    single_mask : bool
+        If set, a single mask will be created, regardless of
+        the number of input integrations. Otherwise, the mask will
+        be a 3D cube, with one plane for each integration.
+    save_mask : bool
+        If set, a mask model is created and returned along with the mask
+        array. If not, the `mask_model` returned is None.
+
+    Returns
+    -------
+    background_mask : array-like of bool
+        Mask array, with True indicating background pixels, False
+        indicating source pixels.
+    mask_model : `~jwst.datamodel.JwstDataModel` or None
+        A datamodel containing the background mask, if `save_mask`
+        is True.
+    """
+    # Check for a user-supplied mask image. If provided, use it.
+    if user_mask is not None:
+        mask_model = datamodels.open(user_mask)
+        background_mask = (mask_model.data.copy()).astype(np.bool_)
+        mask_model.close()
+        del mask_model
+    else:
+        # Create the pixel mask that will be used to indicate which pixels
+        # to include in the 1/f noise measurements. Basically, we're setting
+        # all illuminated pixels to False, so that they do not get used, and
+        # setting all unilluminated pixels to True (so they DO get used).
+        log.info("Creating mask")
+        background_mask = create_mask(
+            image_model,
+            mask_science_regions=mask_science_regions,
+            n_sigma=n_sigma, fit_histogram=fit_histogram,
+            single_mask=single_mask)
+
+    # Store the mask image in a model, if requested
+    if save_mask:
+        mask_model = _make_intermediate_model(image_model, background_mask)
+    else:
+        mask_model = None
+
+    return background_mask, mask_model
+
+
+def _check_data_shapes(input_model, background_mask):
+    """
+    Check data shape for input model and background mask.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodel.JwstDataModel`
+        The input data model.
+    background_mask : array-like of bool
+        The background mask.
+
+    Returns
+    -------
+    mismatch : bool
+        If True, the background mask does not match the data
+        and the step should be skipped.
+    ndim : int
+        Number of dimensions in the input data.
+    nints : int
+        Number of integrations in the input data.
+    ngroups : int
+        Number of groups in the input data.
+    """
+    mismatch = True
+    image_shape = input_model.data.shape[-2:]
+    ndim = input_model.data.ndim
+    if ndim == 2:
+        nints = 1
+        ngroups = 1
+    else:
+        nints = input_model.data.shape[0]
+        if ndim == 3:
+            ngroups = 1
+        else:
+            ngroups = input_model.data.shape[1] - 1
+
+        # Check for 3D mask
+        if background_mask.ndim == 2:
+            if nints > 1:
+                log.info("Data has multiple integrations, but mask is 2D: "
+                         "the same mask will be used for all integrations.")
+        elif background_mask.shape[0] != nints:
+            log.warning("Mask does not match data shape. Step will be skipped.")
+            return mismatch, ndim, nints, ngroups
+
+    # Check that mask matches 2D data shape
+    if background_mask.shape[-2:] != image_shape:
+        log.warning("Mask does not match data shape. Step will be skipped.")
+        return mismatch, ndim, nints, ngroups
+
+    return False, ndim, nints, ngroups
+
+
+def _clean_one_image(image, mask, background_method, background_box_size,
+                     n_sigma, fit_method, detector, fc, axis_to_correct,
+                     fit_by_channel):
+    """
+    Clean an image by fitting and removing background noise.
+
+    Parameters
+    ----------
+    image : array-like of float
+        The input image.
+    mask : array-like of bool
+        The scene mask.  All non-background signal should be
+        marked as True.
+    background_method : str
+        The method for fitting the background.
+    background_box_size : tuple of int
+        Background box size, used with `background_method`
+        is 'model'.
+    n_sigma : float
+        N-sigma rejection level for outliers.
+    fit_method : str
+        The method for fitting the noise after background
+        signal is removed.
+    detector : str
+        The detector name.
+    fc : tuple of int
+        Frequency cutoff values, used for `fit_method` = 'fft'.
+    axis_to_correct : int
+        Axis along which noise appears, used for
+        `fit_method` = 'median'.
+    fit_by_channel : bool
+        Option to fit noise in detector channels separately,
+        used for `fit_method` = 'median'
+
+    Returns
+    -------
+    cleaned_image : array-like of float
+        The cleaned image.
+    background : array-like of float
+        The background fit and removed prior to cleaning.
+        Used for diagnostic purposes.
+    success : bool
+        True if cleaning proceeded as expected; False if
+        cleaning failed and the step should be skipped.
+    """
+    success = True
+
+    # Make sure data is float32
+    image = image.astype(np.float32)
+
+    # Find and replace/mask NaNs
+    nan_pix = np.isnan(image)
+    image[nan_pix] = 0.0
+    mask[nan_pix] = False
+
+    # If there's no good data remaining in the mask,
+    # skip this image
+    if not np.any(mask):
+        return None, None, success
+
+    # Fit and remove a background level
+    if str(background_method).lower() == 'none':
+        background = 0.0
+        bkg_sub = image
+    else:
+        background = background_level(
+            image, mask, background_method=background_method,
+            background_box_size=background_box_size)
+        log.debug(f'Background level: {np.nanmedian(background):.5g}')
+        bkg_sub = image - background
+
+        # Flag more signal in the background subtracted image,
+        # with sigma set by the lower half of the distribution only
+        clip_to_background(
+            bkg_sub, mask, sigma_lower=n_sigma,
+            sigma_upper=n_sigma, lower_half_only=True)
+
+    # Clean the noise
+    if fit_method == 'fft':
+        if bkg_sub.shape == (2048, 2048):
+            # Clean a full-frame image
+            cleaned_image = fft_clean_full_frame(bkg_sub, mask, detector)
+        else:
+            # Clean a subarray image
+            cleaned_image = fft_clean_subarray(bkg_sub, mask, detector, fc=fc)
+    else:
+        cleaned_image = median_clean(bkg_sub, mask, axis_to_correct,
+                                     fit_by_channel=fit_by_channel)
+    if cleaned_image is None:
+        return None, None, False
+
+    # Restore the background level
+    cleaned_image += background
+
+    # Restore NaNs in cleaned image
+    cleaned_image[nan_pix] = np.nan
+
+    return cleaned_image, background, success
+
+
+def _mask_unusable(mask, dq):
+    """
+    Mask unusable data, according to DQ.
+
+    Currently, JUMP and DO_NOT_USE flags are used.
+    The mask is updated in place
+
+    Parameters
+    ----------
+    mask : array-like of bool
+        Input mask, updated in place.
+    dq : array-like of int
+        DQ flag array matching the mask shape.
+    """
+    dnu = (dq & dqflags.group['DO_NOT_USE']) > 0
+    mask[dnu] = False
+    jump = (dq & dqflags.group['JUMP_DET']) > 0
+    mask[jump] = False
+
+
+def do_correction(input_model, input_dir=None, fit_method='median',
+                  fit_by_channel=False, background_method='median',
+                  background_box_size=None,
+                  mask_science_regions=False, n_sigma=2.0, fit_histogram=False,
+                  single_mask=True, user_mask=None, save_mask=False,
+                  save_background=False, save_noise=False):
+    """
+    Apply the 1/f noise correction.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodel.JwstDataModel`
+        Science data to be corrected.
+
+    input_dir : str
+        Path to the input directory.  Used by sub-steps (e.g. assign_wcs
+        for NIRSpec MOS data) to find auxiliary data.
+
+    fit_method : {'fft', 'median'}, optional
+        The algorithm to use to fit background noise.
+
+    fit_by_channel : bool, optional
+        If set, flicker noise is fit independently for each detector channel.
+        Ignored for MIRI, for subarray data, and for `fit_method` = 'fft'.
+
+    background_method : {'median', 'model', None}, optional
+        If 'median', the preliminary background to remove and restore
+        is a simple median of the background data.  If 'model', the
+        background data is fit with a low-resolution model via
+        `~photutils.background.Background2D`.  If None, the background
+        value is 0.0.
+
+    background_box_size : tuple of int, optional
+        Box size for the data grid used by `Background2D` when
+        `background_method` = 'model'. For best results, use a box size
+        that evenly divides the input image shape.
+
+    mask_science_regions : bool, optional
+        For NIRSpec, mask regions of the image defined by WCS bounding
+        boxes for slits/slices, as well as any regions known to be
+        affected by failed-open MSA shutters.  For MIRI imaging, mask
+        regions of the detector not used for science.
+
+    n_sigma : float, optional
+        N-sigma rejection level for finding outliers.
+
+    fit_histogram : bool, optional
+        If set, the 'sigma' used with `n_sigma` for clipping outliers
+        is derived from a Gaussian fit to a histogram of values.
+        Otherwise, a simple iterative sigma clipping is performed.
+
+    single_mask : bool, optional
+        If set, a single mask will be created, regardless of
+        the number of input integrations. Otherwise, the mask will
+        be a 3D cube, with one plane for each integration.
+
+    user_mask : str or None, optional
+        Path to user-supplied mask image.
+
+    save_mask : bool, optional
+        Switch to indicate whether the mask should be saved.
+
+    save_background : bool, optional
+        Switch to indicate whether the fit background should be saved.
+
+    save_noise : bool, optional
+        Switch to indicate whether the fit noise should be saved.
+
+    Returns
+    -------
+    output_model : `~jwst.datamodel.JwstDataModel`
+        Corrected data.
+
+    mask_model : `~jwst.datamodel.JwstDataModel`
+        Pixel mask to be saved or None.
+
+    background_model : `~jwst.datamodel.JwstDataModel`
+        Background model to be saved or None.
+
+    noise_model : `~jwst.datamodel.JwstDataModel`
+        Background model to be saved or None.
+
+    status : {'COMPLETE', 'SKIPPED'}
+        Completion status.  If errors were encountered, status = 'SKIPPED'
+        and the output data matches the input data.  Otherwise,
+        status = 'COMPLETE'.
+    """
+    # Track the completion status, for various failure conditions
+    status = 'SKIPPED'
+
+    detector = input_model.meta.instrument.detector.upper()
+    subarray = input_model.meta.subarray.name.upper()
+    exp_type = input_model.meta.exposure.type
+    slowaxis = input_model.meta.subarray.slowaxis
+    log.info(f'Input exposure type is {exp_type}, detector={detector}')
+
+    # Check for a valid input that we can work on
+    if not _check_input(exp_type, fit_method):
+        return input_model, None, None, None, status
+
+    output_model = input_model.copy()
+
+    # Get parameters needed for subsequent corrections, as appropriate
+    # to the input data
+    axis_to_correct, background_method, fit_by_channel, fc = _standardize_parameters(
+            exp_type, subarray, slowaxis, background_method, fit_by_channel)
+
+    # Make a rate file if needed
+    if user_mask is None:
+        image_model = _make_processed_rate_image(
+            input_model, single_mask, input_dir, exp_type, mask_science_regions)
+    else:
+        image_model = input_model
+
+    # Make a mask model from the user input or the rate data
+    background_mask, mask_model = _make_scene_mask(
+        user_mask, image_model, mask_science_regions,
+        n_sigma, fit_histogram, single_mask, save_mask)
+
+    log.info(f"Cleaning image {input_model.meta.filename}")
+
+    # Check data shapes for 2D, 3D, or 4D inputs
+    mismatch, ndim, nints, ngroups = _check_data_shapes(
+        input_model, background_mask)
+    if mismatch:
+        return output_model, None, None, None, status
+
+    # Close the draft rate model if created - it is no longer needed.
+    if image_model is not input_model:
+        image_model.close()
+        del image_model
+
+    # Make a background cube for saving, if desired
+    if save_background:
+        background_to_save = np.zeros_like(input_model.data)
+    else:
+        background_to_save = None
+
+    # Loop over integrations and groups (even if there's only 1)
+    for i in range(nints):
+        log.debug(f"Working on integration {i + 1}")
+        for j in range(ngroups):
+            log.debug(f"Working on group {j + 1}")
+
+            # Copy the scene mask, for further flagging
+            if background_mask.ndim == 3:
+                mask = background_mask[i].copy()
+            else:
+                mask = background_mask.copy()
+
+            # Get the relevant image data
+            if ndim == 2:
+                image = input_model.data
+            elif ndim == 3:
+                image = input_model.data[i]
+            else:
+                # Ramp data input:
+                # subtract the current group from the next one
+                image = input_model.data[i, j+1] - input_model.data[i, j]
+                dq = input_model.groupdq[i, j+1]
+
+                # Mask any DNU and JUMP pixels
+                _mask_unusable(mask, dq)
+
+            # Clean the image
+            cleaned_image, background, success = _clean_one_image(
+                image, mask, background_method, background_box_size, n_sigma,
+                fit_method, detector, fc, axis_to_correct, fit_by_channel)
+
+            if not success:
+                # Cleaning failed for internal reasons - probably the
+                # mask is not a good match to the data.
+                log.error(f'Cleaning failed for integration {i + 1}, group {j + 1}')
+
+                # Restore input data to make sure any partial changes
+                # are thrown away
+                output_model.data = input_model.data.copy()
+                return output_model, None, None, None, status
+
+            if cleaned_image is None:
+                # Cleaning did not proceed because the image is bad:
+                # leave it as is but continue correcting the rest
+                log.warning(f'No usable data in integration {i+1}, group {j+1}. '
+                            f'Skipping correction for this image.')
+                continue
+
+            # Store the cleaned image in the output model
+            if ndim == 2:
+                output_model.data = cleaned_image
+                if save_background:
+                    background_to_save[:] = background
+            elif ndim == 3:
+                output_model.data[i] = cleaned_image
+                if save_background:
+                    background_to_save[i] = background
+            else:
+                # Add the cleaned data diff to the previously cleaned group,
+                # rather than the noisy input group
+                output_model.data[i, j+1] = output_model.data[i, j] + cleaned_image
+                if save_background:
+                    background_to_save[i, j+1] = background
+
+    # Store the background image in a model, if requested
+    if save_background:
+        background_model = _make_intermediate_model(output_model, background_to_save)
+    else:
+        background_model = None
+
+    # Make a fit noise model for diagnostic purposes by
+    # diffing the input and output models
+    if save_noise:
+        noise_data = output_model.data - input_model.data
+        noise_model = _make_intermediate_model(output_model, noise_data)
+    else:
+        noise_model = None
+
+    # Set completion status
+    status = 'COMPLETE'
+
+    return output_model, mask_model, background_model, noise_model, status
diff --git a/jwst/clean_flicker_noise/clean_flicker_noise_step.py b/jwst/clean_flicker_noise/clean_flicker_noise_step.py
new file mode 100644
index 000000000..a61c179bc
--- /dev/null
+++ b/jwst/clean_flicker_noise/clean_flicker_noise_step.py
@@ -0,0 +1,146 @@
+from stdatamodels.jwst import datamodels
+
+from ..stpipe import Step
+from . import clean_flicker_noise
+
+__all__ = ["CleanFlickerNoiseStep"]
+
+
+class CleanFlickerNoiseStep(Step):
+    """
+    CleanFlickerNoiseStep: This step performs flicker noise correction ("cleaning").
+
+    Input data is expected to be a ramp file (RampModel), in between
+    jump and ramp fitting steps, or a rate file (ImageModel or CubeModel).
+
+    Correction algorithms implemented are:
+        - `fft`: Background noise is fit in frequency space.
+           Implementation is based on the NSClean algorithm, developed
+           by Bernard Rauscher.
+        - `median`: Background noise is characterized by a median
+          along the detector slow axis. Implementation is based on the
+          `image1overf` algorithm, developed by Chris Willott.
+    """
+
+    class_alias = "clean_flicker_noise"
+
+    spec = """
+        fit_method = option('fft', 'median', default='median')  # Noise fitting algorithm
+        fit_by_channel = boolean(default=False)  # Fit noise separately by amplifier (NIR only)
+        background_method = option('median', 'model', None, default='median') # Background fitting algorithm
+        background_box_size = int_list(min=2, max=2, default=None)  # Background box size for modeled background
+        mask_science_regions = boolean(default=False)  # Mask known science regions
+        n_sigma = float(default=2.0)  # Clipping level for non-background signal
+        fit_histogram = boolean(default=False)  # Fit a value histogram to derive sigma
+        single_mask = boolean(default=True)  # Make a single mask for all integrations
+        user_mask = string(default=None)  # Path to user-supplied mask
+        save_mask = boolean(default=False)  # Save the created mask
+        save_background = boolean(default=False)  # Save the fit background
+        save_noise = boolean(default=False)  # Save the fit noise
+        skip = boolean(default=True)  # By default, skip the step
+    """
+
+    def process(self, input):
+        """
+        Fit and subtract 1/f background noise from a ramp data set.
+
+        Parameters
+        ----------
+        input : DataModel
+            Input datamodel to be corrected
+
+        fit_method : str, optional
+            The noise fitting algorithm to use.  Options are 'fft' and 'median'.
+
+        fit_by_channel : bool, optional
+            If set, flicker noise is fit independently for each detector channel.
+            Ignored for MIRI, for subarray data, and for `fit_method` = 'fft'
+
+        background_method : {'median', 'model', None}
+            If 'median', the preliminary background to remove and restore
+            is a simple median of the background data.  If 'model', the
+            background data is fit with a low-resolution model via
+            `~photutils.background.Background2D`.  If None, the background
+            value is 0.0.
+
+        background_box_size : tuple of int, optional
+            Box size for the data grid used by `Background2D` when
+            `background_method` = 'model'. For best results, use a box size
+            that evenly divides the input image shape.
+
+        mask_science_regions : bool, optional
+            For NIRSpec, mask regions of the image defined by WCS bounding
+            boxes for slits/slices, as well as any regions known to be
+            affected by failed-open MSA shutters.  For MIRI imaging, mask
+            regions of the detector not used for science.
+
+        n_sigma : float, optional
+            Sigma clipping threshold to be used in detecting outliers in the image.
+
+        fit_histogram : bool, optional
+            If set, the 'sigma' used with `n_sigma` for clipping outliers
+            is derived from a Gaussian fit to a histogram of values.
+            Otherwise, a simple iterative sigma clipping is performed.
+
+        single_mask : bool, optional
+            If set, a single mask will be created, regardless of
+            the number of input integrations. Otherwise, the mask will
+            be a 3D cube, with one plane for each integration.
+
+        user_mask : None, str, or `~jwst.datamodels.ImageModel`
+            Optional user-supplied mask image; path to file or opened datamodel.
+
+        save_mask : bool, optional
+            Save the computed mask image.
+
+        save_background : bool, optional
+            Save the computed background image.
+
+        save_noise : bool, optional
+            Save the computed noise image.
+
+        Returns
+        -------
+        output_model : DataModel
+            The flicker noise corrected datamodel
+        """
+
+        # Open the input data model
+        with datamodels.open(input) as input_model:
+
+            result = clean_flicker_noise.do_correction(
+                input_model, self.input_dir, self.fit_method, self.fit_by_channel,
+                self.background_method, self.background_box_size,
+                self.mask_science_regions, self.n_sigma, self.fit_histogram,
+                self.single_mask, self.user_mask,
+                self.save_mask, self.save_background, self.save_noise)
+            output_model, mask_model, background_model, noise_model, status = result
+
+            # Save the mask, if requested
+            if self.save_mask and mask_model is not None:
+                mask_path = self.make_output_path(
+                    basepath=input_model.meta.filename, suffix='mask')
+                self.log.info(f"Saving mask file {mask_path}")
+                mask_model.save(mask_path)
+                mask_model.close()
+
+            # Save the background, if requested
+            if self.save_background and background_model is not None:
+                bg_path = self.make_output_path(
+                    basepath=input_model.meta.filename, suffix='flicker_bkg')
+                self.log.info(f"Saving background file {bg_path}")
+                background_model.save(bg_path)
+                background_model.close()
+
+            # Save the noise, if requested
+            if self.save_noise and noise_model is not None:
+                noise_path = self.make_output_path(
+                    basepath=input_model.meta.filename, suffix='flicker_noise')
+                self.log.info(f"Saving noise file {noise_path}")
+                noise_model.save(noise_path)
+                noise_model.close()
+
+            # Set the step completion status
+            output_model.meta.cal_step.clean_flicker_noise = status
+
+        return output_model
diff --git a/jwst/nsclean/lib.py b/jwst/clean_flicker_noise/lib.py
similarity index 99%
rename from jwst/nsclean/lib.py
rename to jwst/clean_flicker_noise/lib.py
index d4bc05c0f..0af311dd2 100644
--- a/jwst/nsclean/lib.py
+++ b/jwst/clean_flicker_noise/lib.py
@@ -517,7 +517,6 @@ class NSCleanSubarray:
         if return_fit:
             return(rfft)
 
-
     def clean(self, weight_fit=True, return_model=False):
         """
         Clean the data
@@ -533,7 +532,7 @@ class NSCleanSubarray:
 
         Returns
         -------
-        data : 2D float array
+        data : array-like of float
             The cleaned data array.
         """ 
         self.fit(weight_fit=weight_fit)  # Fit the background model
diff --git a/jwst/lib/suffix.py b/jwst/lib/suffix.py
index 3fff82ba9..a7af86beb 100644
--- a/jwst/lib/suffix.py
+++ b/jwst/lib/suffix.py
@@ -159,6 +159,8 @@ _calculated_suffixes = {
  'spec3pipeline',
  'nsclean',
  'nscleanstep',
+ 'clean_flicker_noise',
+ 'cleanflickernoisestep',
  'outlierdetectionstep',
  'groupscalestep',
  'hlspstep',
diff --git a/jwst/msaflagopen/msaflag_open.py b/jwst/msaflagopen/msaflag_open.py
index b348ff4fb..39a49e7b7 100644
--- a/jwst/msaflagopen/msaflag_open.py
+++ b/jwst/msaflagopen/msaflag_open.py
@@ -84,10 +84,10 @@ def flag(input_datamodel, failed_slitlets, wcs_refnames):
         science data with DQ flags of affected modified
 
     """
-    #
     # Use the machinery in assign_wcs to create a WCS object for the bad shutters
     pipeline = slitlets_wcs(input_datamodel, wcs_refnames, failed_slitlets)
     wcs = WCS(pipeline)
+
     # Create output as a copy of the input science data model so we can overwrite
     # the wcs with the wcs for the failed open shutters
     # Have to make sure the EXP_TYPE is NRS_MSASPEC so that nrs_wcs_set_input works,
@@ -98,31 +98,30 @@ def flag(input_datamodel, failed_slitlets, wcs_refnames):
     temporary_copy.meta.exposure.type = 'NRS_MSASPEC'
     dq_array = input_datamodel.dq
     for slitlet in failed_slitlets:
-        #
         # Pick the WCS for this slitlet from the WCS of the exposure
         thiswcs = nrs_wcs_set_input(temporary_copy, slitlet.name)
-        #
+
         # Convert the bounding box for this slitlet to a set of indices to use as a slice
         xmin, xmax, ymin, ymax = boundingbox_to_indices(temporary_copy,
                                                         thiswcs.bounding_box)
-        #
+
         # Make a grid of points within the slice
         y_indices, x_indices = np.mgrid[ymin:ymax, xmin:xmax]
-        #
+
         # Calculate the arrays of coordinates for each pixel in the slice
         coordinate_array = thiswcs(x_indices, y_indices)
-        #
+
         # The coordinate_array is a tuple of arrays, one for each output coordinate
         # In this case there should be 3 arrays, one each for RA, Dec and Wavelength
         # For pixels outside the slitlet, the arrays have NaN
-        #
+
         # Make a subarray from these coordinate arrays by setting pixels that aren't
         # NaN to FAILDOPENFLAG, the rest to 0
         dq_subarray = wcs_to_dq(coordinate_array, FAILEDOPENFLAG)
-        #
+
         # bitwise-or this subarray with the slice in the original exposure's DQ array
         dq_array = or_subarray_with_array(dq_array, dq_subarray, xmin, xmax, ymin, ymax)
-    #
+
     # Set the dq array of the input datamodel to the corrected dq array
     input_datamodel.dq = dq_array
     return input_datamodel
@@ -148,7 +147,7 @@ def boundingbox_to_indices(data_model, bounding_box):
         Range of indices of overlap between science data array and bounding box
 
     """
-    nrows, ncols = data_model.data.shape
+    nrows, ncols = data_model.data.shape[-2:]
     ((x1, x2), (y1, y2)) = bounding_box
     xmin = int(min(x1, x2))
     xmin = max(xmin, 0)
@@ -243,5 +242,5 @@ def or_subarray_with_array(dq_array, dq_subarray, xmin, xmax, ymin, ymax):
     Bitwise-or the slice of the dq array with the section corresponding to the
     failed open shutter
     """
-    dq_array[ymin:ymax, xmin:xmax] = np.bitwise_or(dq_array[ymin:ymax, xmin:xmax], dq_subarray)
+    dq_array[..., ymin:ymax, xmin:xmax] = np.bitwise_or(dq_array[..., ymin:ymax, xmin:xmax], dq_subarray)
     return dq_array
diff --git a/jwst/nsclean/nsclean.py b/jwst/nsclean/nsclean.py
deleted file mode 100644
index 27c650954..000000000
--- a/jwst/nsclean/nsclean.py
+++ /dev/null
@@ -1,508 +0,0 @@
-import logging
-import numpy as np
-import gwcs
-from astropy.stats import sigma_clipped_stats
-
-from jwst import datamodels
-from jwst.assign_wcs import nirspec
-from stdatamodels.jwst.datamodels import dqflags
-
-from jwst.nsclean.lib import NSClean, NSCleanSubarray
-
-log = logging.getLogger(__name__)
-log.setLevel(logging.DEBUG)
-
-
-def mask_ifu_slices(input_model, mask):
-
-    """Find pixels located within IFU slices and flag them in the
-    mask, so that they do not get used.
-
-    Parameters
-    ----------
-    input_model : data model object
-        science data
-
-    mask : 2D bool array
-        input mask that will be updated
-
-    Returns
-    -------
-    mask : 2D bool array
-        output mask with additional flags for science pixels
-    """
-
-    log.info("Finding slice pixels for an IFU image")
-
-    # Initialize global DQ map to all zero (OK to use)
-    dqmap = np.zeros_like(input_model.dq)
-
-    # Get the wcs objects for all IFU slices
-    list_of_wcs = nirspec.nrs_ifu_wcs(input_model)
-
-    # Loop over the IFU slices, finding the valid region for each
-    for (k, ifu_wcs) in enumerate(list_of_wcs):
-
-        # Construct array indexes for pixels in this slice
-        x, y = gwcs.wcstools.grid_from_bounding_box(ifu_wcs.bounding_box,
-                                                    step=(1, 1),
-                                                    center=True)
-        # Get the world coords for all pixels in this slice;
-        # all we actually need are wavelengths
-        coords = ifu_wcs(x, y)
-        dq = dqmap[y.astype(int), x.astype(int)]
-        wl = coords[2]
-        # set non-NaN wavelength locations as do not use (one)
-        valid = ~np.isnan(wl)
-        dq = dq[valid]
-        x = x[valid]
-        y = y[valid]
-        dq[:] = 1
-
-        # Copy DQ for this slice into global DQ map
-        dqmap[y.astype(int), x.astype(int)] = dq
-
-    # Now set all non-zero locations in the mask to False (do not use)
-    mask[dqmap==1] = False
-
-    return mask
-
-
-def mask_slits(input_model, mask):
-
-    """Find pixels located within MOS or fixed slit footprints
-    and flag them in the mask, so that they do not get used.
-
-    Parameters
-    ----------
-    input_model : data model object
-        science data
-
-    mask : 2D bool array
-        input mask that will be updated
-
-    Returns
-    -------
-    mask : 2D bool array
-        output mask with additional flags for slit pixels
-    """
-
-    from jwst.extract_2d.nirspec import offset_wcs
-
-    log.info("Finding slit/slitlet pixels")
-
-    # Get the slit-to-msa frame transform from the WCS object
-    slit2msa = input_model.meta.wcs.get_transform('slit_frame', 'msa_frame')
-
-    # Loop over the slits, marking all the pixels within each bounding
-    # box as False (do not use) in the mask.
-    # Note that for 3D masks (TSO mode), all planes will be set to the same value.
-    for slit in slit2msa.slits:
-        slit_wcs = nirspec.nrs_wcs_set_input(input_model, slit.name)
-        xlo, xhi, ylo, yhi = offset_wcs(slit_wcs)
-        mask[..., ylo:yhi, xlo:xhi] = False
-
-    return mask
-
-
-def create_mask(input_model, mask_spectral_regions, n_sigma):
-    """Create the pixel mask needed for setting which pixels to use
-    for measuring 1/f noise.
-
-    Parameters
-    ----------
-    input_model : data model object
-        science data
-
-    mask_spectral_regions : bool
-        mask slit/slice regions defined in WCS
-
-    n_sigma : float
-        sigma threshold for masking outliers
-
-    Returns
-    -------
-    mask : 2D or 3D bool array
-        image mask
-
-    nan_pix : array
-        indexes of image locations with NaN values
-    """
-    exptype = input_model.meta.exposure.type.lower()
-
-    # Initialize mask to all True. Subsequent operations will mask
-    # out pixels that contain signal.
-    # Note: mask will be 3D for BOTS mode data
-    mask = np.full(np.shape(input_model.dq), True)
-
-    # Mask any reference pixels
-    ref_pix = input_model.dq & dqflags.pixel['REFERENCE_PIXEL']
-    mask[ref_pix > 0] = False
-    
-    # If IFU, mask all pixels contained in the IFU slices
-    if exptype == 'nrs_ifu' and mask_spectral_regions:
-        mask = mask_ifu_slices(input_model, mask)
-
-    # If MOS or FS, mask all pixels affected by open slitlets
-    if exptype in ['nrs_fixedslit', 'nrs_brightobj', 'nrs_msaspec'] and mask_spectral_regions:
-        mask = mask_slits(input_model, mask)
-
-    # If IFU or MOS, mask pixels affected by failed-open shutters
-    if exptype in ['nrs_ifu', 'nrs_msaspec']:
-        open_pix = input_model.dq & dqflags.pixel['MSA_FAILED_OPEN']
-        mask[open_pix > 0] = False
-
-    # Temporarily reset NaN pixels and mask them.
-    # Save the list of NaN pixel coords, so that they can be reset at the end.
-    nan_pix = np.isnan(input_model.data)
-    input_model.data[nan_pix] = 0
-    mask[nan_pix] = False
-
-    # If IFU or MOS, mask the fixed-slit area of the image; uses hardwired indexes
-    if exptype == 'nrs_ifu':
-        log.info("Masking the fixed slit region for IFU data.")
-        mask[922:1116, :] = False
-    elif exptype == 'nrs_msaspec':
-        # check for any slits defined in the fixed slit quadrant:
-        # if there is nothing there of interest, mask the whole FS region
-        slit2msa = input_model.meta.wcs.get_transform('slit_frame', 'msa_frame')
-        is_fs = [s.quadrant == 5 for s in slit2msa.slits]
-        if not any(is_fs):
-            log.info("Masking the fixed slit region for MOS data.")
-            mask[922:1116, :] = False
-        else:
-            log.info("Fixed slits found in MSA definition; "
-                     "not masking the fixed slit region for MOS data.")
-
-    # Mask outliers using sigma clipping stats.
-    # For BOTS mode, which uses 3D data, loop over each integration separately.
-    if len(input_model.data.shape) == 3:
-        for i in range(input_model.data.shape[0]):
-            _, median, sigma = sigma_clipped_stats(input_model.data[i], mask=~mask[i], mask_value=0, sigma=5.0)
-            outliers = input_model.data[i] > (median + n_sigma * sigma)
-            mask[i][outliers] = False
-    else:
-        _, median, sigma = sigma_clipped_stats(input_model.data, mask=~mask, mask_value=0, sigma=5.0)
-        outliers = input_model.data > (median + n_sigma * sigma)
-        mask[outliers] = False
-
-
-    # Return the mask and the record of which pixels were NaN in the input;
-    # it'll be needed later
-    return mask, nan_pix
-
-
-def clean_full_frame(detector, image, mask):
-    """Clean a full-frame (2048x2048) image.
-
-    Parameters
-    ----------
-    detector : str
-        The name of the detector from which the data originate.
-
-    image : 2D float array
-        The image to be cleaned.
-
-    mask : 2D bool array
-        The mask that indicates which pixels are to be used in fitting.
-
-    Returns
-    -------
-    cleaned_image : 2D float array
-        The cleaned image.
-    """
-
-    # Instantiate the cleaner
-    cleaner = NSClean(detector, mask)
-
-    # Clean the image
-    try:
-        cleaned_image = cleaner.clean(image, buff=True)
-    except np.linalg.LinAlgError:
-        log.warning("Error cleaning image; step will be skipped")
-        return None
-
-    return cleaned_image
-
-
-def clean_subarray(detector, image, mask, npix_iter=512,
-                   fc=(1061, 1211, 49943, 49957),
-                   exclude_outliers=True, sigrej=4, minfrac=0.05):
-    """Clean a subarray image.
-
-    Parameters
-    ----------
-    detector : str
-        The name of the detector from which the data originate.
-
-    image : 2D float array
-        The image to be cleaned.
-
-    mask : 2D bool array
-        The mask that indicates which pixels are to be used in fitting.
-
-    npix_iter : int
-        Number of pixels to process simultaneously.  Default 512.  Should
-        be at least a few hundred to access sub-kHz frequencies in areas
-        where most pixels are available for fitting.  Previous default
-        behavior corresponds to npix_iter of infinity.
-
-    fc : tuple
-        Apodizing filter definition. These parameters are tunable. The
-        defaults happen to work well for NIRSpec BOTS exposures.
-          1) Unity gain for f < fc[0]
-          2) Cosine roll-off from fc[0] to fc[1]
-          3) Zero gain from fc[1] to fc[2]
-          4) Cosine roll-on from fc[2] to fc[3]
-        Default (1061, 1211, 49943, 49957)
-    
-    exclude_outliers : bool
-        Find and mask outliers in the fit?  Default True
-
-    sigrej : float
-        Number of sigma to clip when identifying outliers.  Default 4.
-
-    minfrac : float
-        Minimum fraction of pixels locally available in the mask in
-        order to attempt a correction.  Default 0.05 (i.e., 5%).
-        
-    Returns
-    -------
-    cleaned_image : 2D float array
-        The cleaned image.
-    """
-
-    # Flip the image to detector coords. NRS1 requires a transpose
-    # of the axes, while NRS2 requires a transpose and flip.
-    if detector == "NRS1":
-        image = image.transpose()
-        mask = mask.transpose()
-    else:
-        image = image.transpose()[::-1]
-        mask = mask.transpose()[::-1]
-
-    # We must do the masking of discrepant pixels here: it just
-    # doesn't work if we wait and do it in the cleaner.  This is
-    # basically copied from lib.py.  Use a robust estimator for
-    # standard deviation, then exclude discrepant pixels and their
-    # four nearest neighbors from the fit.
-    
-    if exclude_outliers:
-        med = np.median(image[mask])
-        std = 1.4825 * np.median(np.abs((image - med)[mask]))
-        outlier = mask & (np.abs(image - med) > sigrej * std)
-        
-        mask = mask&(~outlier)
-        
-        # also get four nearest neighbors of flagged pixels
-        mask[1:] = mask[1:] & (~outlier[:-1])
-        mask[:-1] = mask[:-1] & (~outlier[1:])
-        mask[:, 1:] = mask[:, 1:] & (~outlier[:, :-1])
-        mask[:, :-1] = mask[:, :-1] & (~outlier[:, 1:])
-
-    # Used to determine the fitting intervals along the slow scan
-    # direction.  Pre-pend a zero so that sum_mask[i] is equal
-    # to np.sum(mask[:i], axis=1).
-
-    sum_mask = np.array([0] + list(np.cumsum(np.sum(mask, axis=1))))
-
-    # i1 will be the first row with a nonzero element in the mask
-    # imax will be the last row with a nonzero element in the mask
-    
-    nonzero_mask_element = np.sum(mask, axis=1) > 0
-    i1 = np.amin(np.arange(mask.shape[0])[nonzero_mask_element])
-    imax = np.amax(np.arange(mask.shape[0])[nonzero_mask_element])
-        
-    i1_vals = []
-    di_list = []
-    models = []
-    while i1 <= imax:
-
-        # Want npix_iter available pixels in this section.  If
-        # there are fewer than 1.5*npix_iter available pixels in
-        # the rest of the image, just go to the end.
-        
-        for k in range(i1 + 1, imax + 2):
-            if (sum_mask[k] - sum_mask[i1] > npix_iter
-                    and sum_mask[-1] - sum_mask[i1] > 1.5 * npix_iter):
-                break
-            
-        di = k - i1
-
-        i1_vals += [i1]
-        di_list += [di]
-
-        # Fit this section only if at least minpct% of the pixels
-        # are available for finding the background.  Don't flag
-        # outliers section-by-section; we have to do that earlier
-        # over the full array to get reliable values for the mean
-        # and standard deviation.
-        
-        if np.mean(mask[i1:i1 + di]) > minfrac:
-            cleaner = NSCleanSubarray(image[i1:i1 + di], mask[i1:i1 + di],
-                                      fc=fc, exclude_outliers=False)
-            models += [cleaner.clean(return_model=True)]
-        else:
-            log.warning("Insufficient reference pixels for NSClean around "
-                        "row %d; no correction will be made here." % i1)
-            models += [np.zeros(image[i1:i1 + di].shape)]
-
-        # If we have reached the end of the array, we are finished.
-        if k == imax + 1:
-            break
-
-        # Step forward by half an interval so that we have
-        # overlapping fitting regions.
-        
-        i1 += max(int(np.round(di/2)), 1)
-            
-    model = np.zeros(image.shape)
-    tot_wgt = np.zeros(image.shape)
-
-    # When we combine different corrections computed over
-    # different intervals, each one the highest weight towards the
-    # center of its interval and less weight towards the edge.
-    # Use nonzero weights everywhere so that if only one
-    # correction is available it gets unit weight when we
-    # normalize.
-    
-    for i in range(len(models)):
-        wgt = 1.001 - np.abs(np.linspace(-1, 1, di_list[i]))[:, np.newaxis]
-        model[i1_vals[i]:i1_vals[i] + di_list[i]] += wgt*models[i]
-        tot_wgt[i1_vals[i]:i1_vals[i] + di_list[i]] += wgt
-
-    # don't divide by zero
-    tot_wgt[model == 0] = 1
-    model /= tot_wgt
-    cleaned_image = image - model
-
-    # Restore the cleaned image to the science frame
-    if detector == "NRS1":
-        cleaned_image = cleaned_image.transpose()
-    else:
-        cleaned_image = cleaned_image[::-1].transpose()
-
-    return cleaned_image
-
-
-def do_correction(input_model, mask_spectral_regions, n_sigma, save_mask, user_mask):
-
-    """Apply the NSClean 1/f noise correction
-
-    Parameters
-    ----------
-    input_model : data model object
-        science data to be corrected
-
-    mask_spectral_regions : bool
-        Mask slit/slice regions defined in WCS
-
-    n_sigma : float
-        n-sigma rejection level for finding outliers
-
-    save_mask : bool
-        switch to indicate whether the mask should be saved
-
-    user_mask : str or None
-        Path to user-supplied mask image
-
-    Returns
-    -------
-    output_model : `~jwst.datamodel.JwstDataModel`
-        corrected data
-
-    mask_model : `~jwst.datamodel.JwstDataModel`
-        pixel mask to be saved or None
-    """
-
-    detector = input_model.meta.instrument.detector.upper()
-    exp_type = input_model.meta.exposure.type
-    log.info(f'Input exposure type is {exp_type}, detector={detector}')
-
-    output_model = input_model.copy()
-
-    # Check for a user-supplied mask image. If so, use it.
-    if user_mask is not None:
-        mask_model = datamodels.open(user_mask)
-        Mask = (mask_model.data.copy()).astype(np.bool_)
-
-        # Reset and save list of NaN pixels in the input image
-        nan_pix = np.isnan(input_model.data)
-        input_model.data[nan_pix] = 0
-        Mask[nan_pix] = False
-
-    else:
-        # Create the pixel mask that'll be used to indicate which pixels
-        # to include in the 1/f noise measurements. Basically, we're setting
-        # all illuminated pixels to False, so that they do not get used, and
-        # setting all unilluminated pixels to True (so they DO get used).
-        # For BOTS mode the mask will be 3D, to accommodate changes in masked
-        # pixels per integration.
-        log.info("Creating mask")
-        Mask, nan_pix = create_mask(input_model, mask_spectral_regions, n_sigma)
-
-        # Store the mask image in a model, if requested
-        if save_mask:
-            if len(Mask.shape) == 3:
-                mask_model = datamodels.CubeModel(data=Mask)
-            else:
-                mask_model = datamodels.ImageModel(data=Mask)
-        else:
-            mask_model = None
-
-    log.info(f"Cleaning image {input_model.meta.filename}")
-
-    # Setup for handling 2D or 3D inputs
-    if len(input_model.data.shape) == 3:
-        nints = input_model.data.shape[0]
-        # Check for 3D mask
-        if len(Mask.shape) == 2:
-            log.warning("Data are 3D, but mask is 2D. Step will be skipped.")
-            output_model.meta.cal_step.nsclean = 'SKIPPED'
-            return output_model, None
-    else:
-        nints = 1
-
-    # Loop over integrations (even if there's only 1)
-    for i in range(nints):
-        log.debug(f" working on integration {i+1}")
-        if len(input_model.data.shape) == 3:
-            image = np.float32(input_model.data[i])
-            mask = Mask[i]
-        else:
-            image = np.float32(input_model.data)
-            mask = Mask
-
-        if input_model.data.shape[-2:] == (2048, 2048):
-            # Clean a full-frame image
-            cleaned_image = clean_full_frame(detector, image, mask)
-        else:
-            # BOTS and ALLSLITS exposures should be fitting different
-            # ranges of 1/f frequencies.  Be less aggressive with
-            # fitting higher frequencies in ALLSLITS mode.
-            if input_model.meta.subarray.name.upper() == "ALLSLITS":
-                fc = (150, 200, 49943, 49957)
-            else:
-                fc = (1061, 1211, 49943, 49957)
-            
-            # Clean a subarray image
-            cleaned_image = clean_subarray(detector, image, mask, fc=fc)
-
-        # Check for failure
-        if cleaned_image is None:
-            output_model.meta.cal_step.nsclean = 'SKIPPED'
-            break
-        else:
-            # Store the cleaned image in the output model
-            if len(output_model.data.shape) == 3:
-                output_model.data[i] = cleaned_image
-            else:
-                output_model.data = cleaned_image
-
-    # Restore NaN's from original image
-    output_model.data[nan_pix] = np.nan
-
-    # Set completion status
-    output_model.meta.cal_step.nsclean = 'COMPLETE'
-
-    return output_model, mask_model
diff --git a/jwst/nsclean/nsclean_step.py b/jwst/nsclean/nsclean_step.py
index 79a61128e..83c25f42f 100644
--- a/jwst/nsclean/nsclean_step.py
+++ b/jwst/nsclean/nsclean_step.py
@@ -1,7 +1,7 @@
 from stdatamodels.jwst import datamodels
 
+from jwst.clean_flicker_noise import clean_flicker_noise
 from ..stpipe import Step
-from . import nsclean
 
 __all__ = ["NSCleanStep"]
 
@@ -10,15 +10,25 @@ class NSCleanStep(Step):
     """
     NSCleanStep: This step performs 1/f noise correction ("cleaning")
     of NIRSpec images, using the "NSClean" method.
+
+    NOTE: This step is a deprecated alias to the ``clean_flicker_noise`` step.
     """
 
     class_alias = "nsclean"
 
     spec = """
-        mask_spectral_regions = boolean(default=True)  # Mask WCS-defined regions
+        fit_method = option('fft', 'median', default='fft')  # Noise fitting algorithm
+        fit_by_channel = boolean(default=False)  # Fit noise separately by amplifier
+        background_method = option('median', 'model', None, default=None) # Background fitting algorithm
+        background_box_size = int_list(min=2, max=2, default=None)  # Background box size for modeled background
+        mask_spectral_regions = boolean(default=True)  # Mask WCS-defined spectral regions
         n_sigma = float(default=5.0)  # Clipping level for outliers
-        save_mask = boolean(default=False)  # Save the created mask
+        fit_histogram = boolean(default=False)  # Fit a value histogram to derive sigma
+        single_mask = boolean(default=False)  # Make a single mask for all integrations
         user_mask = string(default=None)  # Path to user-supplied mask
+        save_mask = boolean(default=False)  # Save the created mask
+        save_background = boolean(default=False)  # Save the fit background
+        save_noise = boolean(default=False)  # Save the fit noise
         skip = boolean(default=True)  # By default, skip the step
     """
 
@@ -29,39 +39,102 @@ class NSCleanStep(Step):
         Parameters
         ----------
         input : `~jwst.datamodels.ImageModel`, `~jwst.datamodels.IFUImageModel`
-            Input datamodel to be corrected
+            Input datamodel to be corrected.
+
+        fit_method : str, optional
+            The background fit algorithm to use.  Options are 'fft' and 'median';
+            'fft' performs the original NSClean implementation.
+
+        fit_by_channel : bool, optional
+            If set, flicker noise is fit independently for each detector channel.
+            Ignored for subarray data and for `fit_method` = 'fft'.
+
+        background_method : {'median', 'model', None}
+            If 'median', the preliminary background to remove and restore
+            is a simple median of the background data.  If 'model', the
+            background data is fit with a low-resolution model via
+            `~photutils.background.Background2D`.  If None, the background
+            value is 0.0.
+
+        background_box_size : tuple of int, optional
+            Box size for the data grid used by `Background2D` when
+            `background_method` = 'model'. For best results, use a box size
+            that evenly divides the input image shape.
+
+        mask_spectral_regions : bool, optional
+            Mask regions of the image defined by WCS bounding boxes for slits/slices.
 
         n_sigma : float, optional
-            Sigma clipping threshold to be used in detecting outliers in the image
+            Sigma clipping threshold to be used in detecting outliers in the image.
 
-        save_mask : bool, optional
-            Save the computed mask image
+        fit_histogram : bool, optional
+            If set, the 'sigma' used with `n_sigma` for clipping outliers
+            is derived from a Gaussian fit to a histogram of values.
+            Otherwise, a simple iterative sigma clipping is performed.
+
+        single_mask : bool, optional
+            If set, a single mask will be created, regardless of
+            the number of input integrations. Otherwise, the mask will
+            be a 3D cube, with one plane for each integration.
 
         user_mask : None, str, or `~jwst.datamodels.ImageModel`
-            Optional user-supplied mask image; path to file or opened datamodel
+            Optional user-supplied mask image; path to file or opened datamodel.
 
-        mask_spectral_regions : bool, optional
-            Mask regions of the image defined by WCS bounding boxes for slits/slices
+        save_mask : bool, optional
+            Save the computed mask image.
+
+        save_background : bool, optional
+            Save the computed background image.
+
+        save_noise : bool, optional
+            Save the computed noise image.
 
         Returns
         -------
         output_model : `~jwst.datamodels.ImageModel`, `~jwst.datamodels.IFUImageModel`
-            The 1/f corrected datamodel
+            The 1/f corrected datamodel.
         """
+        message = ("The 'nsclean' step is a deprecated alias to 'clean_flicker_noise' "
+                   "and will be removed in future builds.")
+        self.log.warning(message)
 
         # Open the input data model
         with datamodels.open(input) as input_model:
 
             # Do the NSClean correction
-            result = nsclean.do_correction(input_model, self.mask_spectral_regions, self.n_sigma,
-                                           self.save_mask, self.user_mask)
-            output_model, mask_model = result
+            result = clean_flicker_noise.do_correction(
+                input_model, self.input_dir, self.fit_method, self.fit_by_channel,
+                self.background_method, self.background_box_size,
+                self.mask_spectral_regions, self.n_sigma, self.fit_histogram,
+                self.single_mask, self.user_mask,
+                self.save_mask, self.save_background, self.save_noise)
+            output_model, mask_model, background_model, noise_model, status = result
 
             # Save the mask, if requested
             if self.save_mask and mask_model is not None:
-                mask_path = self.make_output_path(basepath=input_model.meta.filename, suffix='mask')
+                mask_path = self.make_output_path(
+                    basepath=input_model.meta.filename, suffix='mask')
                 self.log.info(f"Saving mask file {mask_path}")
                 mask_model.save(mask_path)
                 mask_model.close()
 
+            # Save the background, if requested
+            if self.save_background and background_model is not None:
+                bg_path = self.make_output_path(
+                    basepath=input_model.meta.filename, suffix='flicker_bkg')
+                self.log.info(f"Saving background file {bg_path}")
+                background_model.save(bg_path)
+                background_model.close()
+
+            # Save the noise, if requested
+            if self.save_noise and noise_model is not None:
+                noise_path = self.make_output_path(
+                    basepath=input_model.meta.filename, suffix='flicker_noise')
+                self.log.info(f"Saving noise file {noise_path}")
+                noise_model.save(noise_path)
+                noise_model.close()
+
+            # Set the step completion status
+            output_model.meta.cal_step.nsclean = status
+
         return output_model
diff --git a/jwst/pipeline/calwebb_detector1.py b/jwst/pipeline/calwebb_detector1.py
index 16087ee73..78699ff12 100644
--- a/jwst/pipeline/calwebb_detector1.py
+++ b/jwst/pipeline/calwebb_detector1.py
@@ -22,6 +22,7 @@ from ..reset import reset_step
 from ..persistence import persistence_step
 from ..charge_migration import charge_migration_step
 from ..jump import jump_step
+from ..clean_flicker_noise import clean_flicker_noise_step
 from ..ramp_fitting import ramp_fit_step
 from ..gain_scale import gain_scale_step
 
@@ -64,6 +65,7 @@ class Detector1Pipeline(Pipeline):
                  'persistence': persistence_step.PersistenceStep,
                  'charge_migration': charge_migration_step.ChargeMigrationStep,
                  'jump': jump_step.JumpStep,
+                 'clean_flicker_noise': clean_flicker_noise_step.CleanFlickerNoiseStep,
                  'ramp_fit': ramp_fit_step.RampFitStep,
                  'gain_scale': gain_scale_step.GainScaleStep,
                  }
@@ -128,6 +130,9 @@ class Detector1Pipeline(Pipeline):
         # apply the jump step
         input = self.jump(input)
 
+        # apply the clean_flicker_noise step
+        input = self.clean_flicker_noise(input)
+
         # save the corrected ramp data, if requested
         if self.save_calibrated_ramp:
             self.save_model(input, 'ramp')
diff --git a/jwst/pipeline/calwebb_spec2.py b/jwst/pipeline/calwebb_spec2.py
index 37eafd089..83b542ca5 100644
--- a/jwst/pipeline/calwebb_spec2.py
+++ b/jwst/pipeline/calwebb_spec2.py
@@ -249,11 +249,11 @@ class Spec2Pipeline(Pipeline):
 
         # Steps whose order is the same for all types of input:
 
-        # Self-calibrate to flag bad/warm pixels, and apply flags 
+        # Self-calibrate to flag bad/warm pixels, and apply flags
         # to both background and science exposures.
         # skipped by default for all modes
         result = self.badpix_selfcal(
-            calibrated, members_by_type['selfcal'], members_by_type['background'], 
+            calibrated, members_by_type['selfcal'], members_by_type['background'],
             )
         if isinstance(result, datamodels.JwstDataModel):
             # if step is skipped, unchanged sci exposure is returned
diff --git a/jwst/pipeline/clean_noise.cfg b/jwst/pipeline/clean_noise.cfg
new file mode 100644
index 000000000..f52822938
--- /dev/null
+++ b/jwst/pipeline/clean_noise.cfg
@@ -0,0 +1,2 @@
+name = "clean_flicker_noise"
+class = "jwst.clean_flicker_noise.CleanFlickerNoiseStep"
diff --git a/jwst/step.py b/jwst/step.py
index 1c9b67b58..fe051a7b6 100644
--- a/jwst/step.py
+++ b/jwst/step.py
@@ -7,6 +7,7 @@ from .background.background_step import BackgroundStep
 from .badpix_selfcal.badpix_selfcal_step import BadpixSelfcalStep
 from .barshadow.barshadow_step import BarShadowStep
 from .charge_migration.charge_migration_step import ChargeMigrationStep
+from .clean_flicker_noise.clean_flicker_noise_step import CleanFlickerNoiseStep
 from .combine_1d.combine_1d_step import Combine1dStep
 from .coron.stack_refs_step import StackRefsStep
 from .coron.align_refs_step import AlignRefsStep
@@ -76,6 +77,7 @@ __all__ = [
     "AlignRefsStep",
     "KlipStep",
     "HlspStep",
+    "CleanFlickerNoiseStep",
     "CubeBuildStep",
     "CubeSkyMatchStep",
     "DarkCurrentStep",
diff --git a/jwst/stpipe/integration.py b/jwst/stpipe/integration.py
index 597a6afe3..7b787b47a 100644
--- a/jwst/stpipe/integration.py
+++ b/jwst/stpipe/integration.py
@@ -44,6 +44,7 @@ def get_steps():
         ("jwst.step.AlignRefsStep", 'align_refs', False),
         ("jwst.step.KlipStep", 'klip', False),
         ("jwst.step.HlspStep", 'hlsp', False),
+        ("jwst.step.CleanFlickerNoiseStep", 'clean_flicker_noise', False),
         ("jwst.step.CubeBuildStep", 'cube_build', False),
         ("jwst.step.CubeSkyMatchStep", 'cube_skymatch', False),
         ("jwst.step.DarkCurrentStep", 'dark_current', False),


Test patch:
diff --git a/jwst/clean_flicker_noise/tests/__init__.py b/jwst/clean_flicker_noise/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/jwst/clean_flicker_noise/tests/test_clean_flicker_noise.py b/jwst/clean_flicker_noise/tests/test_clean_flicker_noise.py
new file mode 100644
index 000000000..7bc237bc9
--- /dev/null
+++ b/jwst/clean_flicker_noise/tests/test_clean_flicker_noise.py
@@ -0,0 +1,844 @@
+import logging
+
+import gwcs
+import numpy as np
+import pytest
+from stdatamodels.jwst import datamodels
+
+from jwst.assign_wcs.tests.test_nirspec import (
+    create_nirspec_ifu_file, create_nirspec_fs_file)
+from jwst.msaflagopen.tests.test_msa_open import make_nirspec_mos_model, get_file_path
+from jwst.clean_flicker_noise import clean_flicker_noise as cfn
+from jwst.tests.helpers import LogWatcher
+
+
+def add_metadata(model, shape):
+    model.meta.instrument.name = 'MIRI'
+    model.meta.instrument.detector = 'MIRIMAGE'
+    model.meta.instrument.filter = 'F480M'
+    model.meta.observation.date = '2015-10-13'
+    model.meta.observation.time = "00:00:00"
+    model.meta.exposure.type = 'MIR_IMAGE'
+    model.meta.exposure.group_time = 1.0
+    model.meta.subarray.name = 'FULL'
+    model.meta.subarray.xstart = 1
+    model.meta.subarray.ystart = 1
+    model.meta.subarray.xsize = shape[3]
+    model.meta.subarray.ysize = shape[2]
+    model.meta.exposure.frame_time = 1.0
+    model.meta.exposure.ngroups = shape[1]
+    model.meta.exposure.group_time = 1.0
+    model.meta.exposure.nints = shape[0]
+    model.meta.exposure.nframes = 1
+    model.meta.exposure.groupgap = 0
+    model.meta.exposure.readpatt = 'FASTR1'
+    model.meta.subarray.slowaxis = 2
+
+
+def make_small_ramp_model(shape=(3, 5, 10, 10)):
+    rampmodel = datamodels.RampModel(shape)
+    add_metadata(rampmodel, shape)
+
+    # Make data with a constant rate
+    for group in range(shape[1]):
+        rampmodel.data[:, group, :, :] = group
+
+    return rampmodel
+
+
+def make_small_rate_model(shape=(3, 5, 10, 10)):
+    ratemodel = datamodels.ImageModel(shape[2:])
+    add_metadata(ratemodel, shape)
+    ratemodel.data[:] = 1.0
+    return ratemodel
+
+
+def make_small_rateints_model(shape=(3, 5, 10, 10)):
+    ratemodel = datamodels.CubeModel((shape[0], shape[2], shape[3]))
+    add_metadata(ratemodel, shape)
+    ratemodel.data[:] = 1.0
+    return ratemodel
+
+
+def make_nirspec_ifu_model(shape=(2048, 2048)):
+    hdul = create_nirspec_ifu_file(grating='PRISM', filter='CLEAR',
+                                   gwa_xtil=0.35986012, gwa_ytil=0.13448857,
+                                   gwa_tilt=37.1)
+    hdul['SCI'].data = np.ones(shape, dtype=float)
+    rate_model = datamodels.IFUImageModel(hdul)
+    hdul.close()
+    return rate_model
+
+
+def make_nirspec_mos_fs_model():
+    mos_model = make_nirspec_mos_model()
+    mos_model.meta.instrument.msa_metadata_file = get_file_path(
+        'msa_fs_configuration.fits')
+    return mos_model
+
+
+def make_nirspec_fs_model():
+    hdul = create_nirspec_fs_file(grating="G140M", filter="F100LP")
+    hdul['SCI'].data = np.ones((2048, 2048), dtype=float)
+    rate_model = datamodels.ImageModel(hdul)
+    hdul.close()
+
+    # add the slow axis
+    rate_model.meta.subarray.slowaxis = 1
+    return rate_model
+
+
+class MockUpdate:
+    def __init__(self):
+        self.seen = False
+
+    def __call__(self, input_model, mask):
+        self.seen = True
+        return mask
+
+
+@pytest.fixture
+def log_watcher(monkeypatch):
+    # Set a log watcher to check for a log message at any level
+    # in the clean_flicker_noise module
+    watcher = LogWatcher('')
+    logger = logging.getLogger('jwst.clean_flicker_noise.clean_flicker_noise')
+    for level in ['debug', 'info', 'warning', 'error']:
+        monkeypatch.setattr(logger, level, watcher)
+    return watcher
+
+
+def test_make_rate(log_watcher):
+    shape = (3, 5, 10, 10)
+    ramp_model = make_small_ramp_model(shape)
+
+    log_watcher.message = 'Creating draft rate file'
+    result = cfn.make_rate(ramp_model, return_cube=True)
+    assert isinstance(result, datamodels.CubeModel)
+    assert result.data.shape == (shape[0], shape[2], shape[3])
+    assert np.all(result.data == 1.0)
+
+    result = cfn.make_rate(ramp_model, return_cube=False)
+    assert isinstance(result, datamodels.ImageModel)
+    assert result.data.shape == shape[2:]
+    assert np.all(result.data == 1.0)
+
+    ramp_model.close()
+    result.close()
+
+    # Check for expected log message
+    log_watcher.assert_seen()
+
+
+def test_postprocess_rate_nirspec(log_watcher):
+    rate_model = make_nirspec_mos_model()
+
+    log_watcher.message = 'Assigning a WCS'
+    result = cfn.post_process_rate(rate_model, assign_wcs=True)
+    assert isinstance(result.meta.wcs, gwcs.WCS)
+    assert np.sum(result.dq & datamodels.dqflags.pixel['MSA_FAILED_OPEN']) == 0
+    log_watcher.assert_seen()
+
+    log_watcher.message = 'Flagging failed-open'
+    result = cfn.post_process_rate(result, msaflagopen=True)
+    assert np.sum(result.dq & datamodels.dqflags.pixel['MSA_FAILED_OPEN']) > 0
+    log_watcher.assert_seen()
+
+    rate_model.close()
+    result.close()
+
+
+def test_postprocess_rate_miri(log_watcher):
+    rate_model = make_small_rate_model()
+    assert np.sum(rate_model.dq & datamodels.dqflags.pixel['DO_NOT_USE']) == 0
+
+    log_watcher.message = 'Retrieving flat DQ'
+    result = cfn.post_process_rate(rate_model, flat_dq=True)
+    log_watcher.assert_seen()
+    assert np.sum(result.dq & datamodels.dqflags.pixel['DO_NOT_USE']) > 0
+    assert np.all(result.data == rate_model.data)
+
+    rate_model.close()
+    result.close()
+
+
+@pytest.mark.slow
+def test_mask_ifu_slices():
+    rate_model = make_nirspec_ifu_model()
+
+    rate_model = cfn.post_process_rate(rate_model, assign_wcs=True)
+    mask = np.full_like(rate_model.data, True)
+
+    # Mark IFU science regions as False: about 10% of the array
+    cfn.mask_ifu_slices(rate_model, mask)
+    assert np.allclose(np.sum(mask), 0.9 * mask.size, atol=0.01 * mask.size)
+
+    rate_model.close()
+
+
+@pytest.mark.parametrize('exptype,blocked',
+                         [('mos', .012), ('mos_fs', .025), ('fs', .05)])
+def test_mask_slits(exptype, blocked):
+    if exptype == 'mos':
+        rate_model = make_nirspec_mos_model()
+    elif exptype == 'mos_fs':
+        rate_model = make_nirspec_mos_fs_model()
+    else:
+        rate_model = make_nirspec_fs_model()
+
+    # Assign a WCS
+    rate_model = cfn.post_process_rate(rate_model, assign_wcs=True)
+
+    # Mark slits as False
+    mask = np.full_like(rate_model.data, True)
+    cfn.mask_slits(rate_model, mask)
+
+    # Check that the fraction of the array blocked is as expected
+    assert np.allclose(np.sum(mask) / mask.size, 1 - blocked, atol=0.001)
+
+    rate_model.close()
+
+
+@pytest.mark.parametrize('fit_histogram', [True, False])
+@pytest.mark.parametrize('lower_half_only', [True, False])
+def test_clip_to_background(log_watcher, fit_histogram, lower_half_only):
+    """Integrity checks for clipping a simple data array."""
+    # Make an array with normal noise
+    shape = (100, 100)
+    rng = np.random.default_rng(seed=123)
+    image = rng.normal(0, 0.01, size=shape)
+    mask = np.full(shape, True)
+
+    # Add a strong high outlier, strong low outlier
+    image[0, 0] += 1
+    image[1, 1] += -1
+
+    # Center found is close to zero, printed when verbose=True
+    log_watcher.message = "center: 0.000"
+    cfn.clip_to_background(
+        image, mask, fit_histogram=fit_histogram,
+        lower_half_only=lower_half_only, verbose=True)
+    log_watcher.assert_seen()
+
+    # All outliers clipped with defaults, as well as a small
+    # percent of the remaining data
+    assert not mask[0, 0]
+    assert not mask[1, 1]
+    assert np.allclose(np.sum(mask) / mask.size, 0.98, atol=0.019)
+
+    # Upper outlier stays with large sigma_upper
+    mask = np.full(shape, True)
+    cfn.clip_to_background(
+        image, mask, fit_histogram=fit_histogram,
+        lower_half_only=lower_half_only, sigma_upper=1000)
+    assert mask[0, 0]
+    assert not mask[1, 1]
+    assert np.allclose(np.sum(mask) / mask.size, 0.98, atol=0.019)
+
+    # Lower outlier stays with large sigma_lower
+    mask = np.full(shape, True)
+    cfn.clip_to_background(
+        image, mask, fit_histogram=fit_histogram,
+        lower_half_only=lower_half_only, sigma_lower=1000)
+    assert not mask[0, 0]
+    assert mask[1, 1]
+    assert np.allclose(np.sum(mask) / mask.size, 0.98, atol=0.019)
+
+
+def test_clip_to_background_fit_fails(log_watcher):
+    shape = (10, 10)
+
+    # histogram failure: all data NaN
+    log_watcher.message = "Histogram failed"
+    image = np.full(shape, np.nan)
+    mask = np.full(shape, True)
+    with pytest.warns(RuntimeWarning):
+        cfn.clip_to_background(image, mask, fit_histogram=True, verbose=True)
+    assert np.all(mask)
+    log_watcher.assert_seen()
+
+    # if mask is all False, warning is avoided, mask is unchanged
+    mask[:] = False
+    cfn.clip_to_background(image, mask, fit_histogram=True, verbose=True)
+    assert not np.all(mask)
+
+    # fit failure: data is not normal
+    log_watcher.message = "Gaussian fit failed"
+    image = np.full(shape, 0.0)
+    mask = np.full(shape, True)
+    image[5:, 5:] = 0.1
+    cfn.clip_to_background(image, mask, fit_histogram=True, verbose=True)
+    assert np.all(mask)
+    log_watcher.assert_seen()
+
+
+@pytest.mark.parametrize('exptype', ['mos', 'mos_fs', 'ifu'])
+def test_create_mask_nirspec(monkeypatch, exptype):
+    # monkeypatch local functions for speed and check that they are called
+    # (actual behavior tested in separate unit tests)
+    mock = MockUpdate()
+    if exptype == 'mos':
+        # NIRSpec MOS data
+        monkeypatch.setattr(cfn, 'mask_slits', mock)
+        rate_model = make_nirspec_mos_model()
+    elif exptype == 'mos_fs':
+        monkeypatch.setattr(cfn, 'mask_slits', mock)
+        rate_model = make_nirspec_mos_fs_model()
+
+        # also assign a WCS so the FS can be detected
+        rate_model = cfn.post_process_rate(rate_model, assign_wcs=True)
+    else:
+        # NIRSpec IFU data
+        monkeypatch.setattr(cfn, 'mask_ifu_slices', mock)
+        rate_model = make_nirspec_ifu_model((2048, 2038))
+
+    rate_model.data[:] = 1.0
+    mask = cfn.create_mask(rate_model)
+    assert mask.shape == rate_model.data.shape
+
+    # Nothing to mask in uniform data
+    assert np.all(mask)
+
+    # Slit function not called if mask_science_regions not set
+    assert not mock.seen
+
+    # Fixed slit region not blocked
+    assert np.all(mask[cfn.NRS_FS_REGION[0]:cfn.NRS_FS_REGION[1]])
+
+    # Call again but block science regions
+    mask = cfn.create_mask(rate_model, mask_science_regions=True)
+    assert mock.seen
+    if exptype == 'mos_fs':
+        # FS region still not blocked - may need correction
+        assert np.all(mask[cfn.NRS_FS_REGION[0]:cfn.NRS_FS_REGION[1]])
+    else:
+        assert not np.any(mask[cfn.NRS_FS_REGION[0]:cfn.NRS_FS_REGION[1]])
+
+
+def test_create_mask_miri():
+    # small MIRI imaging data
+    rate_model = make_small_rate_model()
+    rate_model.dq[5, 5] = datamodels.dqflags.pixel['DO_NOT_USE']
+
+    mask = cfn.create_mask(rate_model)
+    assert mask.shape == rate_model.data.shape
+
+    # Nothing to mask in uniform data
+    assert np.all(mask)
+
+    # Call again but block non-science regions
+    mask = cfn.create_mask(rate_model, mask_science_regions=True)
+    assert not mask[5, 5]
+    assert np.sum(mask) == mask.size - 1
+
+    rate_model.close()
+
+
+def test_create_mask_from_rateints():
+    # small rateints data
+    rate_model = make_small_rateints_model()
+
+    # Add an outlier in each integration
+    for i in range(rate_model.data.shape[0]):
+        rate_model.data[i, i, i] += 100
+
+    mask = cfn.create_mask(rate_model)
+    assert mask.shape == rate_model.data.shape
+    assert mask.ndim == 3
+
+    # Outlier is masked in each integration
+    assert not mask[0, 0, 0]
+    assert not mask[1, 1, 1]
+    assert not mask[2, 2, 2]
+    assert np.sum(mask) == mask.size - 3
+
+    # Call again but make a single mask:
+    # bad data mask is or-ed across integrations
+    mask = cfn.create_mask(rate_model, single_mask=True)
+    assert mask.shape == rate_model.data.shape[1:]
+    assert mask.ndim == 2
+
+    assert not mask[0, 0]
+    assert not mask[1, 1]
+    assert not mask[2, 2]
+    assert np.sum(mask) == mask.size - 3
+
+    rate_model.close()
+
+
+def test_background_level(log_watcher):
+    shape = (100, 100)
+    image = np.full(shape, 1.0)
+    mask = np.full(shape, True)
+
+    # add an outlier to be clipped
+    image[50, 50] = 1000
+
+    # no background
+    background = cfn.background_level(image, mask, background_method=None)
+    assert background == 0.0
+
+    # median method
+    background = cfn.background_level(image, mask, background_method='median')
+    assert background == 1.0
+
+    # model method
+    background = cfn.background_level(
+        image, mask, background_method='model', background_box_size=(10, 10))
+    assert background.shape == shape
+    assert np.all(background == 1.0)
+
+    # model method with mismatched box size:
+    # warns, but completes successfully
+    log_watcher.message = "does not divide evenly"
+    background = cfn.background_level(
+        image, mask, background_method='model', background_box_size=None)
+    assert background.shape == shape
+    assert np.all(background == 1.0)
+    log_watcher.assert_seen()
+
+    # make image mostly bad, one good region
+    image[:] = np.nan
+    image[20:25, 20:25] = 1.0
+
+    # background fit fails: falls back on simple median
+    log_watcher.message = "Background fit failed, using median"
+    background = cfn.background_level(
+        image, mask, background_method='model', background_box_size=(10, 10))
+    assert background == 1.0
+    log_watcher.assert_seen()
+
+
+@pytest.mark.parametrize('array_type,fraction_good', [('full', 0.65), ('subarray', 0.37)])
+@pytest.mark.parametrize('detector', ['NRS1', 'NRS2'])
+def test_fft_clean(array_type, fraction_good, detector):
+    if array_type == 'full':
+        clean_function = cfn.fft_clean_full_frame
+    else:
+        clean_function = cfn.fft_clean_subarray
+
+    shape = (80, 80)
+    mask = np.full(shape, True)
+
+    # zero image should still come out zero
+    image = np.full(shape, 0.0)
+    cleaned_image = clean_function(image.copy(), mask, detector)
+    assert np.allclose(cleaned_image, 0.0)
+
+    # image with regular vertical pattern, centered on 0
+    high = np.full((80, 16), 0.1)
+    low = np.full((80, 16), -0.1)
+    image = np.hstack([high, low, high, low, high])
+    cleaned_image = clean_function(image.copy(), mask, detector)
+
+    # results should be mostly close to zero,
+    # some artifacts at pattern and edge boundaries,
+    # edge boundaries are currently worse for subarray correction
+    good_correction = np.abs(cleaned_image) < 0.01
+    assert np.allclose(np.sum(good_correction) / cleaned_image.size,
+                       fraction_good, atol=0.1)
+
+
+def test_fft_clean_error(monkeypatch, log_watcher):
+    def raise_error(*args, **kwargs):
+        raise np.linalg.LinAlgError('Linear algebra error')
+
+    clean_function = cfn.fft_clean_full_frame
+    monkeypatch.setattr(cfn.NSClean, 'clean', raise_error)
+
+    shape = (10, 10)
+    mask = np.full(shape, True)
+    image = np.full(shape, 1.0)
+
+    log_watcher.message = "Error cleaning image"
+    cleaned_image = clean_function(image.copy(), mask, 'NRS1')
+    assert cleaned_image is None
+    log_watcher.assert_seen()
+
+
+def test_fft_subarray_clean_error(monkeypatch, log_watcher):
+    shape = (10, 10)
+    image = np.full(shape, 1.0)
+
+    # Mask is all bad: error message, returns None
+    mask = np.full(shape, False)
+    log_watcher.message = "No good pixels"
+    cleaned_image = cfn.fft_clean_subarray(image.copy(), mask, 'NRS1')
+    assert cleaned_image is None
+    log_watcher.assert_seen()
+
+    # Mask is nearly all bad except a pixel at the edge:
+    # triggers a linear algebra error
+    mask[0, 0] = True
+    log_watcher.message = "Error cleaning image"
+    cleaned_image = cfn.fft_clean_subarray(image.copy(), mask, 'NRS1')
+    assert cleaned_image is None
+    log_watcher.assert_seen()
+
+    # Mask is mostly bad: warns but continues
+    mask[5, 5] = True
+    log_watcher.message = "Insufficient reference pixels"
+    cleaned_image = cfn.fft_clean_subarray(image.copy(), mask, 'NRS1', minfrac=0.5)
+    assert np.allclose(cleaned_image, image)
+    log_watcher.assert_seen()
+
+
+def test_median_clean():
+    shape = (20, 20)
+    mask = np.full(shape, True)
+
+    # zero image should still come out zero
+    image = np.full(shape, 0.0)
+    cleaned_image = cfn.median_clean(image, mask, 1)
+    assert np.allclose(cleaned_image, 0.0)
+
+    # image with regular vertical pattern, centered on 0
+    high = np.full((20, 4), 0.1)
+    low = np.full((20, 4), -0.1)
+    image = np.hstack([high, low, high, low, high])
+
+    # clean along vertical axis -
+    # image should be all zero
+    cleaned_image = cfn.median_clean(image, mask, 1)
+    assert np.allclose(cleaned_image, 0.0)
+
+    # clean along horizontal axis -
+    # cleaning removes median value, from high stripes
+    cleaned_image = cfn.median_clean(image, mask, 2)
+    assert np.allclose(cleaned_image, image - 0.1)
+
+
+def test_median_clean_by_channel():
+    shape = (2048, 2048)
+    mask = np.full(shape, True)
+
+    # zero image should still come out zero
+    image = np.full(shape, 0.0)
+    cleaned_image = cfn.median_clean(image, mask, 1, fit_by_channel=True)
+    assert np.allclose(cleaned_image, 0.0)
+
+    # image with regular vertical pattern, centered on 0
+    high = np.full((2048, 16), 0.1)
+    low = np.full((2048, 16), -0.1)
+    image = np.hstack([high, low] * 64)
+
+    # add offset by vertical channel
+    offset = np.zeros_like(image)
+    offset[:512] += 0.2
+    offset[512:1024] += 0.3
+    offset[1024:1536] += 0.4
+    offset[1536:] += 0.5
+    image += offset
+
+    # clean along vertical axis by channel -
+    # image should be all zero
+    cleaned_image = cfn.median_clean(image, mask, 1, fit_by_channel=True)
+    assert np.allclose(cleaned_image, 0.0)
+
+    # clean along vertical axis for the whole image -
+    # vertical pattern is removed, offset by channel stays,
+    # but median value is subtracted
+    cleaned_image = cfn.median_clean(image, mask, 1, fit_by_channel=False)
+    assert np.allclose(cleaned_image, offset - 0.35)
+
+    # clean along horizontal axis -
+    # cleaning removes offset value, leaves vertical pattern
+    cleaned_image = cfn.median_clean(image, mask, 2, fit_by_channel=True)
+    assert np.allclose(cleaned_image, image - offset)
+
+
+@pytest.mark.parametrize('mask_science', [True, False])
+def test_do_correction_miri_imaging_ramp(mask_science):
+    ramp_model = make_small_ramp_model()
+    cleaned, mask, bg, noise, status = cfn.do_correction(
+        ramp_model, mask_science_regions=mask_science)
+
+    # uniform data, correction has no effect
+    assert cleaned.data.shape == ramp_model.shape
+    assert np.allclose(cleaned.data, ramp_model.data)
+    assert status == 'COMPLETE'
+
+    # extra models are not created
+    assert mask is None
+    assert bg is None
+    assert noise is None
+
+    ramp_model.close()
+    cleaned.close()
+
+
+@pytest.mark.parametrize('mask_science', [True, False])
+def test_do_correction_nirspec_rate(mask_science):
+    rate_model = make_nirspec_fs_model()
+    cleaned, mask, bg, noise, status = cfn.do_correction(
+        rate_model, mask_science_regions=mask_science)
+
+    # uniform data, correction has no effect
+    assert cleaned.data.shape == rate_model.data.shape
+    assert np.allclose(cleaned.data, rate_model.data)
+    assert status == 'COMPLETE'
+
+    # extra models are not created
+    assert mask is None
+    assert bg is None
+    assert noise is None
+
+    rate_model.close()
+    cleaned.close()
+
+
+@pytest.mark.parametrize('single_mask', [True, False])
+def test_do_correction_rateints(single_mask):
+    rate_model = make_small_rateints_model()
+    cleaned, mask, bg, noise, status = cfn.do_correction(
+        rate_model, single_mask=single_mask, save_mask=True)
+
+    # uniform data, correction has no effect
+    assert cleaned.data.shape == rate_model.shape
+    assert np.allclose(cleaned.data, rate_model.data)
+    assert status == 'COMPLETE'
+
+    if single_mask:
+        assert mask.shape == cleaned.data.shape[-2:]
+    else:
+        assert mask.shape == cleaned.data.shape
+
+    rate_model.close()
+    cleaned.close()
+    mask.close()
+
+
+def test_do_correction_unsupported(log_watcher):
+    ramp_model = make_small_ramp_model()
+    ramp_model.meta.exposure.type = 'MIR_MRS'
+
+    log_watcher.message = "not supported"
+    cleaned, _, _, _, status = cfn.do_correction(ramp_model)
+    assert cleaned is ramp_model
+    assert status == 'SKIPPED'
+    log_watcher.assert_seen()
+
+    ramp_model.close()
+
+
+def test_do_correction_no_fit_by_channel(log_watcher):
+    ramp_model = make_small_ramp_model()
+
+    # fit_by_channel is only used for NIR data -
+    # log a warning, but step still completes
+    log_watcher.message = "can only be used for full-frame NIR"
+    cleaned, _, _, _, status = cfn.do_correction(ramp_model, fit_by_channel=True)
+    assert status == 'COMPLETE'
+    log_watcher.assert_seen()
+
+    ramp_model.close()
+    cleaned.close()
+
+
+@pytest.mark.parametrize('exptype', ['MIR_IMAGE', 'NRC_IMAGE', 'NIS_IMAGE'])
+def test_do_correction_fft_not_allowed(log_watcher, exptype):
+    ramp_model = make_small_ramp_model()
+    ramp_model.meta.exposure.type = exptype
+
+    # not allowed for MIRI, NIRCAM, NIRISS
+    log_watcher.message = f"cannot be applied to exp_type {exptype}"
+    cleaned, _, _, _, status = cfn.do_correction(ramp_model, fit_method='fft')
+    assert cleaned is ramp_model
+    assert status == 'SKIPPED'
+    log_watcher.assert_seen()
+
+    ramp_model.close()
+
+
+@pytest.mark.parametrize('none_value', [None, 'none', 'None'])
+def test_do_correction_no_background(none_value):
+    ramp_model = make_small_ramp_model()
+
+    # Don't remove background before fitting noise
+    cleaned, _, _, _, _ = cfn.do_correction(
+        ramp_model, background_method=none_value)
+
+    # Output data is all zero: uniform level is removed
+    assert np.allclose(cleaned.data, 0.0)
+
+    ramp_model.close()
+    cleaned.close()
+
+
+@pytest.mark.parametrize('ndim', [2, 3])
+def test_do_correction_user_mask(tmp_path, ndim):
+    ramp_model = make_small_ramp_model()
+
+    if ndim == 3:
+        mask = np.full((ramp_model.shape[0],
+                        ramp_model.shape[2],
+                        ramp_model.shape[3]), True)
+        mask_model = datamodels.CubeModel(mask)
+    else:
+        mask = np.full(ramp_model.shape[-2:], True)
+        mask_model = datamodels.ImageModel(mask)
+    user_mask = str(tmp_path / 'mask.fits')
+    mask_model.save(user_mask)
+    mask_model.close()
+
+    cleaned, output_mask, _, _, _ = cfn.do_correction(
+        ramp_model, user_mask=user_mask, save_mask=True)
+
+    assert np.all(output_mask.data == mask_model.data)
+    assert np.allclose(cleaned.data, ramp_model.data)
+
+    ramp_model.close()
+    cleaned.close()
+    output_mask.close()
+
+
+@pytest.mark.parametrize('input_type', ['rate', 'rateints', 'ramp'])
+def test_do_correction_user_mask_mismatch(tmp_path, input_type, log_watcher):
+    shape = (3, 5, 20, 20)
+    if input_type == 'rate':
+        model = make_small_rate_model(shape)
+    elif input_type == 'rateints':
+        model = make_small_rateints_model(shape)
+    else:
+        model = make_small_ramp_model(shape)
+
+    mask = np.full((10, 10), True)
+    mask_model = datamodels.ImageModel(mask)
+    user_mask = str(tmp_path / 'mask.fits')
+    mask_model.save(user_mask)
+    mask_model.close()
+
+    log_watcher.message = 'Mask does not match'
+    cleaned, output_mask, _, _, status = cfn.do_correction(
+        model, user_mask=user_mask, save_mask=True)
+
+    log_watcher.assert_seen()
+    assert status == 'SKIPPED'
+    assert output_mask is None
+
+    model.close()
+    cleaned.close()
+
+
+def test_do_correction_user_mask_mismatch_integ(tmp_path, log_watcher):
+    shape = (3, 5, 20, 20)
+    model = make_small_rateints_model(shape)
+
+    mask = np.full((2, 20, 20), True)
+    mask_model = datamodels.CubeModel(mask)
+    user_mask = str(tmp_path / 'mask.fits')
+    mask_model.save(user_mask)
+    mask_model.close()
+
+    log_watcher.message = 'Mask does not match'
+    cleaned, output_mask, _, _, status = cfn.do_correction(
+        model, user_mask=user_mask, save_mask=True)
+
+    log_watcher.assert_seen()
+    assert status == 'SKIPPED'
+    assert output_mask is None
+
+    model.close()
+    cleaned.close()
+
+
+@pytest.mark.parametrize('subarray', ['SUBS200A1', 'ALLSLITS'])
+def test_do_correction_fft_subarray(subarray):
+    ramp_model = make_small_ramp_model()
+    ramp_model.meta.exposure.type = 'NRS_FIXEDSLIT'
+    ramp_model.meta.subarray.slowaxis = 1
+    ramp_model.meta.subarray.name = subarray
+
+    cleaned, _, _, noise, _ = cfn.do_correction(
+        ramp_model, fit_method='fft', save_noise=True)
+
+    # uniform data, correction has no effect
+    assert cleaned.data.shape == ramp_model.data.shape
+    assert np.allclose(cleaned.data, ramp_model.data)
+    assert np.allclose(noise.data, 0.0)
+
+    ramp_model.close()
+    cleaned.close()
+    noise.close()
+
+
+def test_do_correction_fft_full():
+    rate_model = make_nirspec_fs_model()
+
+    cleaned, _, _, noise, _ = cfn.do_correction(
+        rate_model, fit_method='fft', save_noise=True)
+
+    # uniform data, correction has no effect
+    assert cleaned.data.shape == rate_model.data.shape
+    assert np.allclose(cleaned.data, rate_model.data)
+    assert np.allclose(noise.data, 0.0)
+
+    rate_model.close()
+    cleaned.close()
+    noise.close()
+
+
+def test_do_correction_clean_fails(monkeypatch, log_watcher):
+    ramp_model = make_small_ramp_model()
+    ramp_model.meta.exposure.type = 'NRS_FIXEDSLIT'
+    ramp_model.meta.subarray.slowaxis = 1
+
+    # Add some noise so that input and output are
+    # expected to be different
+    rng = np.random.default_rng(seed=123)
+    ramp_model.data += rng.normal(0, 0.1, size=ramp_model.data.shape)
+    cleaned, _, _, _, status = cfn.do_correction(
+        ramp_model, fit_method='fft')
+    assert not np.allclose(cleaned.data, ramp_model.data)
+    assert status == 'COMPLETE'
+
+    # Mock a None-value returned from cleaning function
+    monkeypatch.setattr(cfn, 'fft_clean_subarray', lambda *args, **kwargs: None)
+
+    # Call again
+    log_watcher.message = "Cleaning failed"
+    cleaned, _, _, _, status = cfn.do_correction(
+        ramp_model, fit_method='fft')
+
+    # Error message issued, status is skipped,
+    # output data is the same as input
+    log_watcher.assert_seen()
+    assert status == 'SKIPPED'
+    assert np.allclose(cleaned.data, ramp_model.data)
+
+    ramp_model.close()
+    cleaned.close()
+
+
+@pytest.mark.parametrize('save_type', ['noise', 'background'])
+@pytest.mark.parametrize('input_type', ['rate', 'rateints', 'ramp'])
+def test_do_correction_save_intermediate(save_type, input_type):
+    shape = (3, 5, 20, 20)
+    if input_type == 'rate':
+        model = make_small_rate_model(shape)
+    elif input_type == 'rateints':
+        model = make_small_rateints_model(shape)
+    else:
+        model = make_small_ramp_model(shape)
+
+    save_bg = (save_type == 'background')
+    save_noise = (save_type == 'noise')
+    cleaned, _, background, noise, status = cfn.do_correction(
+        model, save_background=save_bg, save_noise=save_noise)
+
+    # Output background model always matches input shape
+    # and datamodel type
+    if save_bg:
+        assert background.data.shape == model.data.shape
+        assert type(background) is type(model)
+        background.close()
+    else:
+        assert background is None
+    if save_noise:
+        assert noise.data.shape == model.data.shape
+        assert type(noise) is type(model)
+        noise.close()
+    else:
+        assert noise is None
+
+    model.close()
diff --git a/jwst/clean_flicker_noise/tests/test_clean_flicker_noise_step.py b/jwst/clean_flicker_noise/tests/test_clean_flicker_noise_step.py
new file mode 100644
index 000000000..d8f2e1dcc
--- /dev/null
+++ b/jwst/clean_flicker_noise/tests/test_clean_flicker_noise_step.py
@@ -0,0 +1,82 @@
+import os
+
+import pytest
+from stdatamodels.jwst import datamodels
+
+from jwst.clean_flicker_noise import CleanFlickerNoiseStep
+from .test_clean_flicker_noise import make_small_ramp_model
+
+
+@pytest.mark.parametrize('skip', [True, False])
+def test_output_type(skip):
+    input_model = make_small_ramp_model()
+    cleaned = CleanFlickerNoiseStep.call(input_model, skip=skip)
+
+    # output is a ramp model either way
+    assert isinstance(cleaned, datamodels.RampModel)
+
+    # status may be SKIPPED or COMPLETE
+    if skip:
+        assert cleaned.meta.cal_step.clean_flicker_noise == 'SKIPPED'
+    else:
+        assert cleaned.meta.cal_step.clean_flicker_noise == 'COMPLETE'
+
+    input_model.close()
+    cleaned.close()
+
+
+def test_save_mask(tmp_path):
+    input_model = make_small_ramp_model()
+    input_model.meta.filename = 'test_jump.fits'
+
+    output_dir = str(tmp_path)
+    cleaned = CleanFlickerNoiseStep.call(
+        input_model, skip=False, save_mask=True, output_dir=output_dir)
+
+    mask_file = str(tmp_path / 'test_mask.fits')
+    assert os.path.isfile(mask_file)
+
+    mask_model = datamodels.open(mask_file)
+    assert isinstance(mask_model, datamodels.ImageModel)
+
+    input_model.close()
+    cleaned.close()
+    mask_model.close()
+
+
+def test_save_background(tmp_path):
+    input_model = make_small_ramp_model()
+    input_model.meta.filename = 'test_jump.fits'
+
+    output_dir = str(tmp_path)
+    cleaned = CleanFlickerNoiseStep.call(
+        input_model, skip=False, save_background=True, output_dir=output_dir)
+
+    bkg_file = str(tmp_path / 'test_flicker_bkg.fits')
+    assert os.path.isfile(bkg_file)
+
+    bkg_model = datamodels.open(bkg_file)
+    assert isinstance(bkg_model, datamodels.RampModel)
+
+    input_model.close()
+    cleaned.close()
+    bkg_model.close()
+
+
+def test_save_noise(tmp_path):
+    input_model = make_small_ramp_model()
+    input_model.meta.filename = 'test_jump.fits'
+
+    output_dir = str(tmp_path)
+    cleaned = CleanFlickerNoiseStep.call(
+        input_model, skip=False, save_noise=True, output_dir=output_dir)
+
+    noise_file = str(tmp_path / 'test_flicker_noise.fits')
+    assert os.path.isfile(noise_file)
+
+    noise_model = datamodels.open(noise_file)
+    assert isinstance(noise_model, datamodels.RampModel)
+
+    input_model.close()
+    cleaned.close()
+    noise_model.close()
diff --git a/jwst/msaflagopen/tests/test_msa_open.py b/jwst/msaflagopen/tests/test_msa_open.py
index 131c5eba4..c94f7af81 100644
--- a/jwst/msaflagopen/tests/test_msa_open.py
+++ b/jwst/msaflagopen/tests/test_msa_open.py
@@ -28,6 +28,55 @@ def get_file_path(filename):
     return os.path.join(data_path, filename)
 
 
+def make_nirspec_mos_model():
+    im = ImageModel((2048, 2048))
+    im.meta.wcsinfo = {
+        'dec_ref': -0.00601415671349804,
+        'ra_ref': -0.02073605215697509,
+        'roll_ref': -0.0,
+        'v2_ref': -453.5134,
+        'v3_ref': -373.4826,
+        'v3yangle': 0.0,
+        'vparity': -1}
+
+    im.meta.instrument = {
+        'detector': 'NRS1',
+        'filter': 'F100LP',
+        'grating': 'G140M',
+        'name': 'NIRSPEC',
+        'gwa_tilt': 37.0610,
+        'gwa_xtilt': 0.0001,
+        'gwa_ytilt': 0.0001,
+        'msa_metadata_id': 12}
+
+    im.meta.observation = {
+        'program_number': '1234',
+        'date': '2016-09-05',
+        'time': '8:59:37'}
+
+    im.meta.exposure = {
+        'duration': 11.805952,
+        'end_time': 58119.85416,
+        'exposure_time': 11.776,
+        'frame_time': 0.11776,
+        'group_time': 0.11776,
+        'groupgap': 0,
+        'integration_time': 11.776,
+        'nframes': 1,
+        'ngroups': 100,
+        'nints': 1,
+        'nresets_between_ints': 0,
+        'nsamples': 1,
+        'readpatt': 'NRSRAPID',
+        'sample_time': 10.0,
+        'start_time': 58119.8333,
+        'type': 'NRS_MSASPEC',
+        'zero_frame': False}
+    im.meta.instrument.msa_metadata_file = get_file_path('msa_configuration.fits')
+    im.meta.dither.position_number = 1
+    return im
+
+
 def test_or_subarray_with_array():
     """test bitwise or with array and subarray."""
 
@@ -137,52 +186,7 @@ def test_boundingbox_from_indices():
 
 
 def test_msaflagopen_step():
-    im = ImageModel((2048, 2048))
-    im.meta.wcsinfo = {
-        'dec_ref': -0.00601415671349804,
-        'ra_ref': -0.02073605215697509,
-        'roll_ref': -0.0,
-        'v2_ref': -453.5134,
-        'v3_ref': -373.4826,
-        'v3yangle': 0.0,
-        'vparity': -1}
-
-    im.meta.instrument = {
-        'detector': 'NRS1',
-        'filter': 'F100LP',
-        'grating': 'G140M',
-        'name': 'NIRSPEC',
-        'gwa_tilt': 37.0610,
-        'gwa_xtilt': 0.0001,
-        'gwa_ytilt': 0.0001,
-        'msa_metadata_id': 12}
-
-    im.meta.observation = {
-        'program_number': '1234',
-        'date': '2016-09-05',
-        'time': '8:59:37'}
-
-    im.meta.exposure = {
-        'duration': 11.805952,
-        'end_time': 58119.85416,
-        'exposure_time': 11.776,
-        'frame_time': 0.11776,
-        'group_time': 0.11776,
-        'groupgap': 0,
-        'integration_time': 11.776,
-        'nframes': 1,
-        'ngroups': 100,
-        'nints': 1,
-        'nresets_between_ints': 0,
-        'nsamples': 1,
-        'readpatt': 'NRSRAPID',
-        'sample_time': 10.0,
-        'start_time': 58119.8333,
-        'type': 'NRS_MSASPEC',
-        'zero_frame': False}
-    im.meta.instrument.msa_metadata_file = get_file_path('msa_configuration.fits')
-    im.meta.dither.position_number = 1
-
+    im = make_nirspec_mos_model()
     im = AssignWcsStep.call(im)
     result = MSAFlagOpenStep.call(im)
 
diff --git a/jwst/pipeline/tests/test_calwebb_spec2.py b/jwst/pipeline/tests/test_calwebb_spec2.py
index d060aa9cd..71a39232f 100644
--- a/jwst/pipeline/tests/test_calwebb_spec2.py
+++ b/jwst/pipeline/tests/test_calwebb_spec2.py
@@ -155,4 +155,4 @@ def test_filenotfounderror_raised(capsys):
 
     # Verify the failure is printed to stderr
     captured = capsys.readouterr()
-    assert 'FileNotFoundError' in captured.err
\ No newline at end of file
+    assert 'FileNotFoundError' in captured.err
diff --git a/jwst/regtest/test_miri_image.py b/jwst/regtest/test_miri_image.py
index 5559daf12..192e6dcec 100644
--- a/jwst/regtest/test_miri_image.py
+++ b/jwst/regtest/test_miri_image.py
@@ -47,6 +47,26 @@ def run_detector1_with_average_dark_current(rtdata_module):
     Step.from_cmdline(args)
 
 
+@pytest.fixture(scope="module")
+def run_detector1_with_clean_flicker_noise(rtdata_module):
+    """Run detector1 pipeline on MIRI imaging data with noise cleaning."""
+    rtdata_module.get_data("miri/image/jw01024002001_02101_00001_mirimage_uncal.fits")
+
+    # Run detector1 pipeline only on one of the _uncal files.
+    # Run optional clean_flicker_noise step,
+    # saving extra outputs and masking to science regions
+    args = ["jwst.pipeline.Detector1Pipeline", rtdata_module.input,
+            "--save_calibrated_ramp=True",
+            "--steps.clean_flicker_noise.skip=False",
+            "--steps.clean_flicker_noise.mask_science_regions=True",
+            "--steps.clean_flicker_noise.save_results=True",
+            "--steps.clean_flicker_noise.save_mask=True",
+            "--steps.clean_flicker_noise.save_background=True",
+            "--steps.clean_flicker_noise.save_noise=True",
+            ]
+    Step.from_cmdline(args)
+
+
 @pytest.fixture(scope="module")
 def run_image2(run_detector1, rtdata_module):
     """Run image2 pipeline on the _rate file, saving intermediate products"""
@@ -114,6 +134,29 @@ def test_miri_image_detector1_with_avg_dark_current(run_detector1_with_average_d
     assert diff.identical, diff.report()
 
 
+@pytest.mark.bigdata
+@pytest.mark.parametrize("suffix",
+                         ["clean_flicker_noise", "mask",
+                          "flicker_bkg", "flicker_noise",
+                          "ramp", "rate", "rateints"])
+def test_miri_image_detector1_with_clean_flicker_noise(
+        run_detector1_with_clean_flicker_noise,
+        rtdata_module, fitsdiff_default_kwargs, suffix):
+    """Test detector1 pipeline for MIRI imaging data with noise cleaning."""
+    rtdata = rtdata_module
+    rtdata.input = "jw01024002001_02101_00001_mirimage_uncal.fits"
+    output = f"jw01024002001_02101_00001_mirimage_{suffix}.fits"
+    rtdata.output = output
+
+    rtdata.get_truth(f"truth/test_miri_image_clean_flicker_noise/{output}")
+
+    # Set tolerances so the file comparisons work across architectures
+    fitsdiff_default_kwargs["rtol"] = 1e-4
+    fitsdiff_default_kwargs["atol"] = 1e-4
+    diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs)
+    assert diff.identical, diff.report()
+
+
 @pytest.mark.bigdata
 @pytest.mark.parametrize("suffix", ["assign_wcs", "flat_field", "cal", "i2d"])
 def test_miri_image_image2(run_image2, rtdata_module, fitsdiff_default_kwargs, suffix):
diff --git a/jwst/regtest/test_nircam_image.py b/jwst/regtest/test_nircam_image.py
index 11283baf0..1de07dbf9 100644
--- a/jwst/regtest/test_nircam_image.py
+++ b/jwst/regtest/test_nircam_image.py
@@ -30,6 +30,24 @@ def run_detector1pipeline(rtdata_module):
     Step.from_cmdline(args)
 
 
+@pytest.fixture(scope="module")
+def run_detector1_with_clean_flicker_noise(rtdata_module):
+    """Run detector1 pipeline on NIRCam imaging data with noise cleaning."""
+    rtdata_module.get_data("nircam/image/jw01538046001_03105_00001_nrcalong_uncal.fits")
+
+    # Run detector1 pipeline only on one of the _uncal files.
+    # Run optional clean_flicker_noise step, saving extra outputs
+    args = ["jwst.pipeline.Detector1Pipeline", rtdata_module.input,
+            "--save_calibrated_ramp=True",
+            "--steps.clean_flicker_noise.skip=False",
+            "--steps.clean_flicker_noise.save_results=True",
+            "--steps.clean_flicker_noise.save_mask=True",
+            "--steps.clean_flicker_noise.save_background=True",
+            "--steps.clean_flicker_noise.save_noise=True",
+            ]
+    Step.from_cmdline(args)
+
+
 @pytest.fixture(scope="module")
 def run_image2pipeline(run_detector1pipeline, rtdata_module):
     """Run calwebb_image2 on NIRCam imaging long data"""
@@ -225,3 +243,26 @@ def test_imaging_distortion(rtdata, fitsdiff_default_kwargs):
     assert_allclose(y, model.meta.wcsinfo.crpix2 - 1)
     assert_allclose(raout, ra)
     assert_allclose(decout, dec)
+
+
+@pytest.mark.bigdata
+@pytest.mark.parametrize("suffix",
+                         ["clean_flicker_noise", "mask",
+                          "flicker_bkg", "flicker_noise",
+                          "ramp", "rate", "rateints"])
+def test_nircam_image_detector1_with_clean_flicker_noise(
+        run_detector1_with_clean_flicker_noise,
+        rtdata_module, fitsdiff_default_kwargs, suffix):
+    """Test detector1 pipeline for NIRCam imaging data with noise cleaning."""
+    rtdata = rtdata_module
+    rtdata.input = "jw01538046001_03105_00001_nrcalong_uncal.fits"
+    output = f"jw01538046001_03105_00001_nrcalong_{suffix}.fits"
+    rtdata.output = output
+
+    rtdata.get_truth(f"truth/test_nircam_image_clean_flicker_noise/{output}")
+
+    # Set tolerances so the file comparisons work across architectures
+    fitsdiff_default_kwargs["rtol"] = 1e-4
+    fitsdiff_default_kwargs["atol"] = 1e-4
+    diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs)
+    assert diff.identical, diff.report()
diff --git a/jwst/regtest/test_niriss_image.py b/jwst/regtest/test_niriss_image.py
index 7e186986a..62b268eb6 100644
--- a/jwst/regtest/test_niriss_image.py
+++ b/jwst/regtest/test_niriss_image.py
@@ -36,14 +36,46 @@ def run_detector1(rtdata_module):
     Step.from_cmdline(args)
 
 
+@pytest.fixture(scope="module")
+def run_detector1_with_clean_flicker_noise(rtdata_module):
+    """Run detector1 pipeline on NIRISS imaging data with noise cleaning."""
+    rtdata_module.get_data("niriss/imaging/jw01094001002_02107_00001_nis_uncal.fits")
+
+    # Run detector1 pipeline only on one of the _uncal files.
+    # Run optional clean_flicker_noise step,
+    # saving extra outputs and masking to science regions
+    args = ["jwst.pipeline.Detector1Pipeline", rtdata_module.input,
+            "--save_calibrated_ramp=True",
+            "--steps.clean_flicker_noise.skip=False",
+            "--steps.clean_flicker_noise.save_results=True",
+            "--steps.clean_flicker_noise.save_mask=True",
+            "--steps.clean_flicker_noise.save_background=True",
+            "--steps.clean_flicker_noise.save_noise=True",
+            ]
+    Step.from_cmdline(args)
+
+
 @pytest.mark.bigdata
 @pytest.mark.parametrize("suffix", ["dq_init", "saturation", "superbias",
                                     "refpix", "linearity", "dark_current",
                                     "charge_migration", "jump", "rate", "rateints"])
 def test_niriss_image_detector1(run_detector1, rtdata_module, fitsdiff_default_kwargs, suffix):
+    """Test detector1 pipeline for NIRISS imaging data with noise cleaning.
+    """
+    truth_dir = 'test_niriss_image'
+    _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix, truth_dir)
+
+
+@pytest.mark.bigdata
+@pytest.mark.parametrize("suffix", ["clean_flicker_noise", "mask",
+                                    "flicker_bkg", "flicker_noise",
+                                    "ramp", "rate", "rateints"])
+def test_niriss_image_detector1_with_clean_flicker_noise(
+        run_detector1_with_clean_flicker_noise, rtdata_module, fitsdiff_default_kwargs, suffix):
     """Regression test of detector1 pipeline performed on NIRISS imaging data.
     """
-    _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix)
+    truth_dir = 'test_niriss_image_clean_flicker_noise'
+    _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix, truth_dir)
 
 
 @pytest.mark.bigdata
@@ -77,14 +109,14 @@ def test_niriss_tweakreg_no_sources(rtdata, fitsdiff_default_kwargs):
             result.shelve(model, modify=False)
 
 
-def _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix):
+def _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix, truth_dir):
     """Assertion helper for the above tests"""
     rtdata = rtdata_module
     rtdata.input = "jw01094001002_02107_00001_nis_uncal.fits"
     output = f"jw01094001002_02107_00001_nis_{suffix}.fits"
     rtdata.output = output
 
-    rtdata.get_truth(f"truth/test_niriss_image/{output}")
+    rtdata.get_truth(f"truth/{truth_dir}/{output}")
 
     # Set tolerances so the crf, rscd and rateints file comparisons work across
     # architectures
diff --git a/jwst/regtest/test_nirspec_irs2_detector1.py b/jwst/regtest/test_nirspec_irs2_detector1.py
index 407a8bd2d..8a8441e26 100644
--- a/jwst/regtest/test_nirspec_irs2_detector1.py
+++ b/jwst/regtest/test_nirspec_irs2_detector1.py
@@ -30,6 +30,26 @@ def run_detector1pipeline(rtdata_module):
     ])
 
 
+@pytest.fixture(scope="module")
+def run_detector1_with_clean_flicker_noise(rtdata_module):
+    """Run detector1 pipeline on NIRSpec IRS2 data with noise cleaning."""
+    rtdata_module.get_data("nirspec/irs2/jw01335004001_03101_00002_nrs2_uncal.fits")
+
+    # Run detector1 pipeline only on one of the _uncal files.
+    # Run optional clean_flicker_noise step,
+    # saving extra outputs and masking to science regions
+    args = ["jwst.pipeline.Detector1Pipeline", rtdata_module.input,
+            "--save_calibrated_ramp=True",
+            "--steps.clean_flicker_noise.skip=False",
+            "--steps.clean_flicker_noise.mask_science_regions=True",
+            "--steps.clean_flicker_noise.save_results=True",
+            "--steps.clean_flicker_noise.save_mask=True",
+            "--steps.clean_flicker_noise.save_background=True",
+            "--steps.clean_flicker_noise.save_noise=True",
+            ]
+    Step.from_cmdline(args)
+
+
 @pytest.mark.bigdata
 @pytest.mark.parametrize("suffix", ['dq_init', 'saturation', 'superbias',
                                     'refpix', 'linearity', 'dark_current', 'jump',
@@ -51,3 +71,24 @@ def test_nirspec_irs2_detector1(run_detector1pipeline, rtdata_module,
     fitsdiff_default_kwargs["atol"] = 1e-4
     diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs)
     assert diff.identical, diff.report()
+
+
+@pytest.mark.bigdata
+@pytest.mark.parametrize("suffix", ["clean_flicker_noise", "mask",
+                                    "flicker_bkg", "flicker_noise",
+                                    "ramp", "rate", "rateints"])
+def test_nirspec_irs2_detector1_with_clean_flicker_noise(
+        run_detector1_with_clean_flicker_noise, rtdata_module,
+        fitsdiff_default_kwargs, suffix):
+    """Test detector1 pipeline for NIRSpec IRS2 data with noise cleaning."""
+    rtdata = rtdata_module
+
+    output_filename = f"jw01335004001_03101_00002_nrs2_{suffix}.fits"
+    rtdata.output = output_filename
+    rtdata.get_truth(f"truth/test_nirspec_irs2_clean_flicker_noise/{output_filename}")
+
+    # Set tolerances so comparisons work across architectures
+    fitsdiff_default_kwargs["rtol"] = 1e-4
+    fitsdiff_default_kwargs["atol"] = 1e-4
+    diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs)
+    assert diff.identical, diff.report()