diff --git a/CHANGES.rst b/CHANGES.rst
index 1e73338259..a78806080e 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -141,6 +141,9 @@ refpix
 - Add option for NIRSpec IRS2 to preserve interleaved reference pixels in the
   output file, for calibration and diagnostic purposes. [#8255]
 
+- Add option to correct for mean reference pixel offsets by amplifier and detector
+  column in NIRSpec IRS2 mode. [#8143]
+
 resample
 --------
 
diff --git a/docs/jwst/refpix/description.rst b/docs/jwst/refpix/description.rst
index 43e657fe81..f5fb2d8701 100644
--- a/docs/jwst/refpix/description.rst
+++ b/docs/jwst/refpix/description.rst
@@ -176,7 +176,7 @@ reference pixel correction for data read out using the IRS2 readout pattern.
 See the JdoxIRS2_ page for for an overview, and see Rauscher2017_ for
 details.
 
-The raw data include both the science data and interspersed reference
+The raw data include both the science data and interleaved reference
 pixel values.  The time to read out the entire detector includes not only
 the time to read each pixel of science ("normal") data and some of the
 reference pixels, but also time for the transition between reading normal
@@ -185,17 +185,14 @@ each row and between frames.  For example, it takes the same length of time
 to jump from reading normal pixels to reading reference pixels as it does
 to read one pixel value, about ten microseconds.
 
-Before subtracting the reference pixel and reference output values from
-the science data, some processing is done on the reference values, and the
-CRDS reference file factors are applied.  IRS2 readout is only used for
-full-frame data, never for subarrays.  The full detector is read out
-by four separate amplifiers simultaneously, and the reference output is
-read at the same time.  Each of these five readouts is the same size,
-640 by 2048 pixels (for IRS2).  If the CRDS reference file includes a
-DQ (data quality) BINTABLE extension, interleaved reference pixel values
-will be set to zero if they are flagged as bad in the DQ extension.
-
-At this point the algorithm looks for intermittently bad (or suspicious)
+IRS2 readout is only used for full-frame data, never for subarrays. The full
+detector is read out by four separate amplifiers simultaneously, and the
+reference output is read at the same time.  Each of these five readouts is
+the same size, 640 by 2048 pixels, each containing a repeating set of 8
+normal pixel readouts, 4 interleaved reference pixel readouts, and 8 more
+normal pixel readouts.
+
+The first step in processing IRS2 data is to look for intermittently bad
 reference pixels. This is done by calculating the means and standard
 deviations per reference pixel column, as well as the absolute value of the
 difference between readout pairs, across all groups within each integration.
@@ -203,7 +200,7 @@ The robust mean and standard deviation of each of these arrays is then
 computed. Values greater than the robust mean plus the standard
 deviation, times a factor to avoid overcorrection, are flagged as bad
 pixels.  Readout pairs are always flagged together, and are flagged across
-all groups and integrations. Bad values are replaced by values from the
+all groups and integrations. Bad values will be replaced by values from the
 nearest reference group within the same amplifier, respecting parity
 (even/oddness).  The replacement value is the average of upper and lower
 values if both are good, or directly using the upper or lower values if only
@@ -212,6 +209,24 @@ good adjacent neighbor that does not match parity, that value is used.  If
 there are no good replacement values, the bad pixel is set to 0.0 to be
 interpolated over in the IRS2 correction to follow.
 
+After flagging bad reference pixels, the step performs an optional
+correction for overall mean reference pixel offsets by amplifier and
+column parity. The algorithm described above for the traditional NIR readout
+mode is applied to IRS2 data to perform this correction, with two small
+differences:
+
+    #. Side pixel correction is never applied for IRS2 data.
+
+    #. "Even" and "odd" refer to detector column addresses, rather than
+       data array locations, to ensure that interleaved reference pixel
+       columns are accounted for correctly.
+
+After the mean offsets are subtracted and bad pixels are replaced, some processing
+is done on the remaining reference values, and the CRDS reference file
+factors are applied. If the CRDS reference file includes a DQ (data quality)
+BINTABLE extension, interleaved reference pixel values will be set to zero if
+they are flagged as bad in the DQ extension.
+
 The next step in this processing is to
 copy the science data and the reference pixel data separately to temporary
 1-D arrays (both of length 712 * 2048); this is done separately for each
diff --git a/jwst/refpix/irs2_subtract_reference.py b/jwst/refpix/irs2_subtract_reference.py
index 5fd33cd0ce..7eaff7627a 100644
--- a/jwst/refpix/irs2_subtract_reference.py
+++ b/jwst/refpix/irs2_subtract_reference.py
@@ -607,37 +607,37 @@ def flag_bad_refpix(datamodel, n_sigma=3.0, flag_only=False, replace_only=False)
                         rp_stds.append(rp_s)
                         rp_diffs.append(rp_d)
 
-            ref_pix = np.array(ref_pix, dtype=int)
-            rp_diffs = np.array(rp_diffs)
-            rp_means = np.array(rp_means)
-            rp_stds = np.array(rp_stds)
-            pair_pixel = ref_pix + 1
-
-            # clipped stats for all tests
-            mean_of_diffs, _, std_of_diffs = sigma_clipped_stats(rp_diffs, sigma=n_sigma)
-            mean_of_means, _, std_of_means = sigma_clipped_stats(rp_means, sigma=n_sigma)
-            mean_of_stds, _, std_of_stds = sigma_clipped_stats(rp_stds, sigma=n_sigma)
-
-            # find the additional intermittent bad pixels, marking both readouts
-            high_diffs = (rp_diffs - mean_of_diffs) > (n_sigma * std_of_diffs)
-            high_means = (rp_means - mean_of_means) > (n_sigma * std_of_means)
-            high_stds = (rp_stds - mean_of_stds) > (n_sigma * std_of_stds)
-
-            log.debug(f'High diffs={np.sum(high_diffs)}, '
-                      f'high means={np.sum(high_means)}, '
-                      f'high stds={np.sum(high_stds)}')
-
-            int_bad[ref_pix[high_diffs]] = True
-            int_bad[pair_pixel[high_diffs]] = True
-            int_bad[ref_pix[high_means]] = True
-            int_bad[pair_pixel[high_means]] = True
-            int_bad[ref_pix[high_stds]] = True
-            int_bad[pair_pixel[high_stds]] = True
-
-            log.debug(f'{np.sum(int_bad[offset:offset + amplifier])} '
-                      f'suspicious bad reference pixels in '
-                      f'amplifier {k}, integration {j}')
-            mask_bad |= int_bad
+            if not replace_only:
+                ref_pix = np.array(ref_pix, dtype=int)
+                rp_diffs = np.array(rp_diffs)
+                rp_means = np.array(rp_means)
+                rp_stds = np.array(rp_stds)
+                pair_pixel = ref_pix + 1
+
+                # clipped stats for all tests
+                mean_of_diffs, _, std_of_diffs = sigma_clipped_stats(rp_diffs, sigma=n_sigma)
+                mean_of_means, _, std_of_means = sigma_clipped_stats(rp_means, sigma=n_sigma)
+                mean_of_stds, _, std_of_stds = sigma_clipped_stats(rp_stds, sigma=n_sigma)
+
+                # find the additional intermittent bad pixels, marking both readouts
+                high_diffs = (rp_diffs - mean_of_diffs) > (n_sigma * std_of_diffs)
+                high_means = (rp_means - mean_of_means) > (n_sigma * std_of_means)
+                high_stds = (rp_stds - mean_of_stds) > (n_sigma * std_of_stds)
+
+                log.debug(f'High diffs={np.sum(high_diffs)}, '
+                          f'high means={np.sum(high_means)}, '
+                          f'high stds={np.sum(high_stds)}')
+                int_bad[ref_pix[high_diffs]] = True
+                int_bad[pair_pixel[high_diffs]] = True
+                int_bad[ref_pix[high_means]] = True
+                int_bad[pair_pixel[high_means]] = True
+                int_bad[ref_pix[high_stds]] = True
+                int_bad[pair_pixel[high_stds]] = True
+
+                log.debug(f'{np.sum(int_bad[offset:offset + amplifier])} '
+                          f'suspicious bad reference pixels in '
+                          f'amplifier {k}, integration {j}')
+                mask_bad |= int_bad
 
         # replace any flagged pixels if desired
         if not flag_only:
@@ -649,7 +649,10 @@ def flag_bad_refpix(datamodel, n_sigma=3.0, flag_only=False, replace_only=False)
                                offset, offset + amplifier,
                                scipix_n, refpix_r)
 
-    log.info(f'Total bad reference pixels: {np.sum(mask_bad)}')
+    if flag_only:
+        log.info(f'Total bad reference pixels flagged: {np.sum(mask_bad)}')
+    else:
+        log.info(f'Total bad reference pixels replaced: {np.sum(mask_bad)}')
     if pixeldq is not None:
         pixeldq[mask_bad] |= (dqflags.pixel['BAD_REF_PIXEL']
                               | dqflags.pixel['DO_NOT_USE'])
diff --git a/jwst/refpix/reference_pixels.py b/jwst/refpix/reference_pixels.py
index b6f2d509a4..734a1e33bd 100644
--- a/jwst/refpix/reference_pixels.py
+++ b/jwst/refpix/reference_pixels.py
@@ -41,13 +41,16 @@
 #
 #  For MIRI subarray exposures, omit the refpix step.
 
+import logging
+from copy import deepcopy
+
 import numpy as np
 from scipy import stats
-import logging
 
 from stdatamodels.jwst.datamodels import dqflags
 
-from ..lib import reffile_utils
+from ..lib import pipe_utils, reffile_utils
+from .irs2_subtract_reference import make_irs2_mask
 
 log = logging.getLogger(__name__)
 log.setLevel(logging.DEBUG)
@@ -75,6 +78,26 @@
                                 'data': (0, 2048, 1536, 2048)}
                           }
 
+# IRS2 sections for NIRSpec have a different size due to the
+# interleaved reference pixels and the reference sector.
+
+IRS2_reference_sections = {'0': {'top': (2044, 2048, 0, 640),
+                                 'bottom': (0, 4, 0, 640),
+                                 'data': (0, 2048, 0, 640)},
+                           'A': {'top': (2044, 2048, 640, 1280),
+                                 'bottom': (0, 4, 640, 1280),
+                                 'data': (0, 2048, 640, 1280)},
+                           'B': {'top': (2044, 2048, 1280, 1920),
+                                 'bottom': (0, 4, 1280, 1920),
+                                 'data': (0, 2048, 1280, 1920)},
+                           'C': {'top': (2044, 2048, 1920, 2560),
+                                 'bottom': (0, 4, 1920, 2560),
+                                 'data': (0, 2048, 1920, 2560)},
+                           'D': {'top': (2044, 2048, 2560, 3200),
+                                 'bottom': (0, 4, 2560, 3200),
+                                 'data': (0, 2048, 2560, 3200)}
+                           }
+
 # Special behavior is requested for NIRSpec subarrays that do not reach
 # detector edges; for these input models, we will assign the top and bottom
 # four rows as reference pixels to better treat pedestal noise issues.
@@ -195,6 +218,9 @@ def __init__(self, input_model,
         self.odd_even_rows = odd_even_rows
         self.bad_reference_pixels = False
 
+        self.reference_sections = None
+        self.amplifiers = 'ABCD'
+
         # Define temp array for processing every group
         self.pixeldq = self.get_pixeldq()
         self.group = None
@@ -415,7 +441,7 @@ def count_good_side_refpixels(self):
         donotuse = dqflags.pixel['DO_NOT_USE']
         ngood = 0
         for amplifier in 'AD':
-            rowstart, rowstop, colstart, colstop = NIR_reference_sections[amplifier]['side']
+            rowstart, rowstop, colstart, colstop = self.reference_sections[amplifier]['side']
             good = np.where(np.bitwise_and(self.pixeldq[rowstart:rowstop, colstart:colstop], donotuse) != donotuse)
             ngood += len(good[0])
         return ngood
@@ -429,8 +455,8 @@ def count_good_top_bottom_refpixels(self):
             log.debug(f"Edgeless subarray {self.subarray} has {ngood} reference pixels.")
         else:
             for edge in ['top', 'bottom']:
-                for amplifier in 'ABCD':
-                    rowstart, rowstop, colstart, colstop = NIR_reference_sections[amplifier][edge]
+                for amplifier in self.amplifiers:
+                    rowstart, rowstop, colstart, colstop = self.reference_sections[amplifier][edge]
                     log.debug(f"Ref sections for {edge} & {amplifier}: {rowstart, rowstop, colstart, colstop}")
                     good = np.where(np.bitwise_and(self.pixeldq[rowstart:rowstop, colstart:colstop], donotuse) != donotuse)
                     ngood += len(good[0])
@@ -443,7 +469,6 @@ class NIRDataset(Dataset):
 
     Parameters
     ----------
-
     input_model: data model object
         Science data model to be corrected
 
@@ -481,19 +506,129 @@ def __init__(self, input_model,
                                          side_gain,
                                          odd_even_rows=False)
 
-#
-#  Even though the recommendation specifies calculating the mean of the
-#  combined top and bottom reference sections, there's a good chance we
-#  might want to calculate them separately
-#
+        # Set appropriate NIR sections
+        self.is_irs2 = pipe_utils.is_irs2(input_model)
+        if self.is_irs2:
+            self.reference_sections = deepcopy(IRS2_reference_sections)
+            self.amplifiers = '0ABCD'
+            self.irs2_odd_mask = self.make_irs2_odd_mask(input_model)
+        else:
+            self.reference_sections = NIR_reference_sections
+            self.irs2_odd_mask = None
+
+    def make_irs2_odd_mask(self, input_model, scipix_n_default=16, refpix_r_default=4):
+        """
+        Make an odd pixel mask for IRS2 mode.
+
+        The even pixel mask can be generated by inverting the odd mask.
+
+        Parameters
+        ----------
+        input_model : DataModel
+            Input model containing data to mask.
+        scipix_n_default : int, optional
+            Number of regular samples before stepping out to collect
+            reference samples.
+        refpix_r_default : int, optional
+            Number of reference samples before stepping back in to collect
+            regular samples.
+
+        Returns
+        -------
+        odd_mask : NDArray
+            Boolean array matching data column size in detector orientation.
+            True identifies all odd pixels (science and reference).
+        """
+        # Get data information from input model.
+        # (y and x here refer to detector orientation, although input
+        # data has not yet been rotated)
+        ny = input_model.data.shape[-1]  # 2048
+        nx = input_model.data.shape[-2]  # 3200
+
+        # Default n=16, r=4
+        scipix_n = input_model.meta.exposure.nrs_normal
+        if scipix_n is None:
+            log.warning("Keyword NRS_NORM not found; using default value %d" %
+                        scipix_n_default)
+            scipix_n = scipix_n_default
+
+        refpix_r = input_model.meta.exposure.nrs_reference
+        if refpix_r is None:
+            log.warning("Keyword NRS_REF not found; using default value %d" %
+                        refpix_r_default)
+            refpix_r = refpix_r_default
+
+        # If these are not set to standard values, the
+        # reference sections values must be changed to match.
+        n_sector = nx // 5
+        areas = ['top', 'bottom', 'data']  # assuming no 'side'
+        if nx != 3200:
+            for i, amplifier in enumerate('0ABCD'):
+                x_start = n_sector * i
+                x_stop = n_sector * (i + 1)
+                for area in areas:
+                    sec = self.reference_sections[amplifier][area]
+                    self.reference_sections[amplifier][area] = (
+                        sec[0], sec[1], x_start, x_stop)
+
+        # Make a column mask that identifies the reference sector and
+        # all interleaved pixels as False, all science pixels
+        # and standard reference pixels as True
+        x_mask = make_irs2_mask(nx, ny, scipix_n, refpix_r)
+
+        # Switch True/False to identify reference pixels instead of
+        # science pixels
+        x_mask = ~x_mask
+
+        # Treat the reference sector like the other sectors
+        x_mask[:n_sector] = x_mask[n_sector: 2 * n_sector]
+
+        # Find even and odd interleaved pixels:
+        # reference pixels come in two pairs, the first set odd
+        # and the second set even, so pixels
+        #    SSSSSSSSrrrrSSSSSSSSSSSSSSSSrrrrSSSS...
+        # have parity:
+        #    --------0011----------------0011----...
+        # for the interleaved pixels, where the traditional pixels
+        # are:
+        #    01010101----0101010101010101----0101...
+
+        # first two pixels are odd
+        even_interleaved = x_mask.copy()
+        even_interleaved[0::4] = False
+        even_interleaved[1::4] = False
+
+        # second two pixels are even
+        odd_interleaved = x_mask.copy()
+        odd_interleaved[2::4] = False
+        odd_interleaved[3::4] = False
+
+        # Make an odd mask for the image columns in detector orientation.
+        # This will be used both for identifying the correct
+        # reference pixels and for applying the correction later.
+        odd_mask = np.full(nx, False)
+        odd_mask[0::2] = True
+        odd_mask[even_interleaved] = False
+        odd_mask[odd_interleaved] = True
+
+        return odd_mask
+
+    #  Even though the recommendation specifies calculating the mean of the
+    #  combined top and bottom reference sections, there's a good chance we
+    #  might want to calculate them separately
     def collect_odd_refpixels(self, group, amplifier, top_or_bottom):
-        """Collect up the reference pixels corresponding to odd-numbered
-        rows (first, third, fifth, etc, corresponding to even array indices)
+        """Collect odd reference pixels.
 
-        Parameters:
-        -----------
+        For traditional readouts, odd pixels correspond to odd-numbered
+        rows (first, third, fifth, etc.), which are even array indices.
 
-        group: NDArray
+        For IRS2 mode, science and traditional reference pixels have
+        the same parity, but interleaved pixels come in two pairs. The
+        first two are odd and the second two are even.
+
+        Parameters
+        ----------
+        group : NDArray
             The group that is being processed
 
         amplifier: string ['A'|'B'|'C'|'D']
@@ -503,9 +638,9 @@ def collect_odd_refpixels(self, group, amplifier, top_or_bottom):
             String corresponding to whether top or bottom reference pixels
             are bing processed
 
-        Returns:
-
-        oddref: NDArray
+        Returns
+        -------
+        oddref : NDArray
             Array containing all the odd reference pixels
 
         odddq: NDArray
@@ -513,20 +648,33 @@ def collect_odd_refpixels(self, group, amplifier, top_or_bottom):
 
         """
         rowstart, rowstop, colstart, colstop = \
-            NIR_reference_sections[amplifier][top_or_bottom]
+            self.reference_sections[amplifier][top_or_bottom]
+
+        # handle interleaved pixels if needed
+        if self.is_irs2:
+            odd_mask = self.irs2_odd_mask[colstart:colstop]
+            oddref = group[rowstart:rowstop, colstart:colstop][:, odd_mask]
+            odddq = self.pixeldq[rowstart:rowstop, colstart:colstop][:, odd_mask]
+        else:
+            oddref = group[rowstart:rowstop, colstart:colstop:2]
+            odddq = self.pixeldq[rowstart:rowstop, colstart:colstop:2]
 
-        oddref = group[rowstart:rowstop, colstart:colstop: 2]
-        odddq = self.pixeldq[rowstart:rowstop, colstart:colstop: 2]
         return oddref, odddq
 
     def collect_even_refpixels(self, group, amplifier, top_or_bottom):
-        """Collect up the reference pixels corresponding to even-numbered
-        rows (second, fourth, sixth, etc, corresponding to odd array indices)
+        """Collect even reference pixels.
 
-        Parameters:
-        -----------
+        For traditional readouts, even pixels correspond to even-numbered
+        rows (second, fourth, sixth, etc.), which are odd array indices.
 
-        group: NDArray
+        For IRS2 mode, science and traditional reference pixels have
+        the same parity, but interleaved pixels come in two pairs. The
+        first two are odd and the second two are even.
+
+        Parameters
+        ----------
+
+        group : NDArray
             The group that is being processed
 
         amplifier: string ['A'|'B'|'C'|'D']
@@ -536,23 +684,28 @@ def collect_even_refpixels(self, group, amplifier, top_or_bottom):
             String corresponding to whether top or bottom reference pixels
             are bing processed
 
-        Returns:
-
-        oddref: NDArray
-            Array containing all the odd reference pixels
-
-        odddq: NDArray
-            Array containing all the odd dq values for those reference pixels
+        Returns
+        -------
+        evenref : NDArray
+            Array containing all the even reference pixels
 
+        evendq : NDArray
+            Array containing all the even dq values for those reference pixels
         """
-
         rowstart, rowstop, colstart, colstop = \
-            NIR_reference_sections[amplifier][top_or_bottom]
-        #
-        # Even columns start on the second column
-        colstart = colstart + 1
-        evenref = group[rowstart:rowstop, colstart:colstop: 2]
-        evendq = self.pixeldq[rowstart:rowstop, colstart:colstop: 2]
+            self.reference_sections[amplifier][top_or_bottom]
+
+        # handle interleaved pixels if needed
+        if self.is_irs2:
+            even_mask = ~self.irs2_odd_mask[colstart:colstop]
+            evenref = group[rowstart:rowstop, colstart:colstop][:, even_mask]
+            evendq = self.pixeldq[rowstart:rowstop, colstart:colstop][:, even_mask]
+        else:
+            # Even columns start on the second column
+            colstart = colstart + 1
+            evenref = group[rowstart:rowstop, colstart:colstop:2]
+            evendq = self.pixeldq[rowstart:rowstop, colstart:colstop:2]
+
         return evenref, evendq
 
     def get_odd_refvalue(self, group, amplifier, top_or_bottom):
@@ -654,7 +807,7 @@ def get_amplifier_refvalue(self, group, amplifier, top_or_bottom):
             return odd, even
         else:
             rowstart, rowstop, colstart, colstop = \
-                NIR_reference_sections[amplifier][top_or_bottom]
+                self.reference_sections[amplifier][top_or_bottom]
             ref = group[rowstart:rowstop, colstart:colstop]
             dq = self.pixeldq[rowstart:rowstop, colstart:colstop]
             mean = self.sigma_clip(ref, dq)
@@ -683,7 +836,7 @@ def get_refvalues(self, group):
         """
 
         refpix = {}
-        for amplifier in 'ABCD':
+        for amplifier in self.amplifiers:
             refpix[amplifier] = {}
             refpix[amplifier]['odd'] = {}
             refpix[amplifier]['even'] = {}
@@ -721,9 +874,9 @@ def do_top_bottom_correction(self, group, refvalues):
         top and bottom reference pixels
 
         """
-        for amplifier in 'ABCD':
+        for amplifier in self.amplifiers:
             datarowstart, datarowstop, datacolstart, datacolstop = \
-                NIR_reference_sections[amplifier]['data']
+                self.reference_sections[amplifier]['data']
             if self.odd_even_columns:
                 oddreftop = refvalues[amplifier]['odd']['top']
                 oddrefbottom = refvalues[amplifier]['odd']['bottom']
@@ -734,12 +887,19 @@ def do_top_bottom_correction(self, group, refvalues):
                 oddrefsignal = self.average_with_None(oddreftop, oddrefbottom)
                 evenrefsignal = self.average_with_None(evenreftop, evenrefbottom)
                 if oddrefsignal is not None and evenrefsignal is not None:
-                    oddslice = (slice(datarowstart, datarowstop, 1),
-                                slice(datacolstart, datacolstop, 2))
-                    evenslice = (slice(datarowstart, datarowstop, 1),
-                                 slice(datacolstart + 1, datacolstop, 2))
-                    group[oddslice] = group[oddslice] - oddrefsignal
-                    group[evenslice] = group[evenslice] - evenrefsignal
+                    if not self.is_irs2:
+                        oddslice = (slice(datarowstart, datarowstop, 1),
+                                    slice(datacolstart, datacolstop, 2))
+                        evenslice = (slice(datarowstart, datarowstop, 1),
+                                     slice(datacolstart + 1, datacolstop, 2))
+                        group[oddslice] = group[oddslice] - oddrefsignal
+                        group[evenslice] = group[evenslice] - evenrefsignal
+                    else:
+                        dataslice = (slice(datarowstart, datarowstop, 1),
+                                     slice(datacolstart, datacolstop, 1))
+                        odd_mask = self.irs2_odd_mask[datacolstart:datacolstop]
+                        group[dataslice][:, odd_mask] -= oddrefsignal
+                        group[dataslice][:, ~odd_mask] -= evenrefsignal
                 else:
                     pass
             else:
@@ -1427,6 +1587,8 @@ def __init__(self, input_model,
                                           side_gain=False,
                                           odd_even_rows=odd_even_rows)
 
+        self.reference_sections = MIR_reference_sections
+
     def DMS_to_detector(self, integration, group):
         #
         # MIRI data doesn't need transforming
@@ -1463,7 +1625,7 @@ def collect_odd_refpixels(self, group, amplifier, left_or_right):
 
         """
 
-        rowstart, rowstop, column = MIR_reference_sections[amplifier][left_or_right]
+        rowstart, rowstop, column = self.reference_sections[amplifier][left_or_right]
         oddref = group[rowstart:rowstop:2, column]
         odddq = self.pixeldq[rowstart:rowstop:2, column]
         return oddref, odddq
@@ -1494,7 +1656,7 @@ def collect_even_refpixels(self, group, amplifier, left_or_right):
 
         """
 
-        rowstart, rowstop, column = MIR_reference_sections[amplifier][left_or_right]
+        rowstart, rowstop, column = self.reference_sections[amplifier][left_or_right]
         #
         # Even reference pixels start on the second row
         rowstart = rowstart + 1
@@ -1604,7 +1766,7 @@ def get_amplifier_refvalue(self, group, amplifier, left_or_right):
                 self.bad_reference_pixels = True
             return odd, even
         else:
-            rowstart, rowstop, column = MIR_reference_sections[amplifier][left_or_right]
+            rowstart, rowstop, column = self.reference_sections[amplifier][left_or_right]
             ref = group[rowstart:rowstop, column]
             dq = self.pixeldq[rowstart:rowstop, column]
             mean = self.sigma_clip(ref, dq)
@@ -1633,7 +1795,7 @@ def get_refvalues(self, group):
         """
 
         refpix = {}
-        for amplifier in 'ABCD':
+        for amplifier in self.amplifiers:
             refpix[amplifier] = {}
             refpix[amplifier]['odd'] = {}
             refpix[amplifier]['even'] = {}
@@ -1672,9 +1834,9 @@ def do_left_right_correction(self, group, refvalues):
 
         """
 
-        for amplifier in 'ABCD':
+        for amplifier in self.amplifiers:
             datarowstart, datarowstop, datacolstart, datacolstop, stride = \
-                MIR_reference_sections[amplifier]['data']
+                self.reference_sections[amplifier]['data']
             if self.odd_even_rows:
                 oddrefleft = refvalues[amplifier]['odd']['left']
                 oddrefright = refvalues[amplifier]['odd']['right']
diff --git a/jwst/refpix/refpix_step.py b/jwst/refpix/refpix_step.py
index 4baf2a21a2..1e1df78328 100644
--- a/jwst/refpix/refpix_step.py
+++ b/jwst/refpix/refpix_step.py
@@ -24,6 +24,7 @@ class RefPixStep(Step):
         odd_even_rows = boolean(default=True) # Compute reference signal separately for even- and odd-numbered rows
         ovr_corr_mitigation_ftr = float(default=3.0) # Factor to avoid overcorrection of bad reference pixels for IRS2
         preserve_irs2_refpix = boolean(default=False) # Preserve reference pixels in output
+        irs2_mean_subtraction = boolean(default=False) # Apply a mean offset subtraction before IRS2 correction
     """
 
     reference_file_types = ['refpix']
@@ -38,7 +39,20 @@ def process(self, input):
                 # Flag bad reference pixels first
                 datamodel = input_model.copy()
                 irs2_subtract_reference.flag_bad_refpix(
-                    datamodel, n_sigma=self.ovr_corr_mitigation_ftr)
+                    datamodel, n_sigma=self.ovr_corr_mitigation_ftr, flag_only=True)
+
+                # If desired, do the normal refpix correction before IRS2, without
+                # side pixel handling
+                if self.irs2_mean_subtraction:
+                    if self.use_side_ref_pixels:
+                        self.log.info('Turning off side pixel correction for IRS2')
+                        self.use_side_ref_pixels = False
+                    reference_pixels.correct_model(
+                        datamodel, self.odd_even_columns, self.use_side_ref_pixels,
+                        self.side_smoothing_length, self.side_gain, self.odd_even_rows)
+
+                # Now that values are updated, replace bad reference pixels
+                irs2_subtract_reference.flag_bad_refpix(datamodel, replace_only=True)
 
                 # Get the necessary refpix reference file for IRS2 correction
                 self.irs2_name = self.get_reference_file(datamodel, 'refpix')
diff --git a/jwst/refpix/tests/test_refpix.py b/jwst/refpix/tests/test_refpix.py
index 086bd9eb39..f7d5fccaeb 100644
--- a/jwst/refpix/tests/test_refpix.py
+++ b/jwst/refpix/tests/test_refpix.py
@@ -203,10 +203,10 @@ def test_odd_even():
 
 
 @pytest.mark.parametrize('detector,ysize,odd_even',
-                         [('NRS1', 2048, True),
-                          ('NRS1', 2048, False),
-                          ('NRS2', 2048, True),
-                          ('NRS2', 2048, False),])
+                         [('NRS1', 2048, True), ('NRS1', 3200, True),
+                          ('NRS1', 2048, False), ('NRS1', 3200, False),
+                          ('NRS2', 2048, True), ('NRS2', 3200, True),
+                          ('NRS2', 2048, False), ('NRS2', 3200, False)])
 def test_odd_even_amp_nirspec(detector, ysize, odd_even):
     """Check that odd/even columns are applied when flag is set"""
 
@@ -249,7 +249,8 @@ def test_odd_even_amp_nirspec(detector, ysize, odd_even):
     im.pixeldq[:, -4:] = dqflags.pixel['REFERENCE_PIXEL']
 
     # run the step
-    out = RefPixStep.call(im, use_side_ref_pixels=False, odd_even_columns=odd_even)
+    out = RefPixStep.call(im, use_side_ref_pixels=False, odd_even_columns=odd_even,
+                          irs2_mean_subtraction=True)
 
     # values should be different by amp and by odd/even row if specified
     # pick a random pixel to test