diff --git a/CHANGES.rst b/CHANGES.rst
index 48d9eda8c..f89ba5b71 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,7 +1,10 @@
 1.5.3 (unreleased)
 ==================
 
-- 
+ramp_fitting
+------------
+- Added the avg_dark_current parameter that is used to include the
+  dark current in the calculation of the Poisson noise [#236]
 
 1.5.2 (2023-12-13)
 ==================
diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py
index d6aeb3adc..d15dda4d1 100644
--- a/src/stcal/jump/jump.py
+++ b/src/stcal/jump/jump.py
@@ -1,11 +1,13 @@
 import logging
 import multiprocessing
 import time
-
-import cv2 as cv
 import numpy as np
-from astropy import stats
-from astropy.convolution import Ring2DKernel, convolve
+import cv2 as cv
+import astropy.stats as stats
+
+from astropy.convolution import Ring2DKernel
+from astropy.convolution import convolve
+
 
 from . import constants
 from . import twopoint_difference as twopt
@@ -848,6 +850,7 @@ def find_faint_extended(
     Total number of showers detected.
 
     """
+    log.info('Flagging Showers')
     read_noise_2 = readnoise_2d**2
     data = indata.copy()
     data[gdq == sat_flag] = np.nan
diff --git a/src/stcal/ramp_fitting/ols_fit.py b/src/stcal/ramp_fitting/ols_fit.py
index f2c6d5a23..34d4dacc7 100644
--- a/src/stcal/ramp_fitting/ols_fit.py
+++ b/src/stcal/ramp_fitting/ols_fit.py
@@ -19,7 +19,8 @@
 log.setLevel(logging.DEBUG)
 
 
-def ols_ramp_fit_multi(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, max_cores):
+def ols_ramp_fit_multi(
+        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, max_cores):
     """
     Setup the inputs to ols_ramp_fit with and without multiprocessing. The
     inputs will be sliced into the number of cores that are being used for
@@ -57,6 +58,9 @@ def ols_ramp_fit_multi(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, wei
         'half', and 'all'. This is the fraction of cores to use for multi-proc. The
         total number of cores includes the SMT cores (Hyper Threading for Intel).
 
+    avg_dark_current: float (electrons)
+        The average dark current for this detector
+
     Returns
     -------
     image_info : tuple
@@ -96,15 +100,16 @@ def ols_ramp_fit_multi(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, wei
     if number_slices == 1:
         # Single threaded computation
         image_info, integ_info, opt_info = ols_ramp_fit_single(
-            ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting
-        )
+            ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting)
         if image_info is None or integ_info is None:
             return None, None, None
 
         return image_info, integ_info, opt_info
 
     # Call ramp fitting for multi-processor (multiple data slices) case
-    image_info, integ_info, opt_info = ols_ramp_fit_multiprocessing(
+
+    else:
+        image_info, integ_info, opt_info = ols_ramp_fit_multiprocessing(
         ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, number_slices
     )
 
@@ -144,6 +149,8 @@ def ols_ramp_fit_multiprocessing(
     number_slices: int
         The number of slices to partition the data into for multiprocessing.
 
+    avg_dark_current: float
+        The average dark current in electrons
     Return
     ------
     image_info: tuple
@@ -157,8 +164,8 @@ def ols_ramp_fit_multiprocessing(
     """
     log.info("Number of processors used for multiprocessing: %s", number_slices)
     slices, rows_per_slice = compute_slices_for_starmap(
-        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, number_slices
-    )
+        ramp_data, buffsize, save_opt,
+        readnoise_2d, gain_2d, weighting, number_slices)
 
     pool = Pool(processes=number_slices)
     pool_results = pool.starmap(ols_ramp_fit_single, slices)
@@ -432,8 +439,8 @@ def get_max_segs_crs(pool_results):
 
 
 def compute_slices_for_starmap(
-    ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, number_slices
-):
+        ramp_data, buffsize, save_opt,
+        readnoise_2d, gain_2d, weighting, avg_dark_current, number_slices):
     """
     Creates the slices needed for each process for multiprocessing.  The slices
     for the arguments needed for ols_ramp_fit_single.
@@ -473,9 +480,11 @@ def compute_slices_for_starmap(
     start_row = 0
     for k in range(len(rslices)):
         ramp_slice = slice_ramp_data(ramp_data, start_row, rslices[k])
-        rnoise_slice = readnoise_2d[start_row : start_row + rslices[k], :].copy()
-        gain_slice = gain_2d[start_row : start_row + rslices[k], :].copy()
-        slices.insert(k, (ramp_slice, buffsize, save_opt, rnoise_slice, gain_slice, weighting))
+
+        rnoise_slice = readnoise_2d[start_row:start_row + rslices[k], :].copy()
+        gain_slice = gain_2d[start_row:start_row + rslices[k], :].copy()
+        slices.insert(
+            k, (ramp_slice, buffsize, save_opt, rnoise_slice, gain_slice, weighting, avg_dark_current))
         start_row = start_row + rslices[k]
 
     return slices, rslices
@@ -618,7 +627,8 @@ def find_0th_one_good_group(ramp_data):
     ramp_data.one_groups_time = (ramp_data.nframes + 1) * ramp_data.frame_time / 2
 
 
-def ols_ramp_fit_single(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting):
+def ols_ramp_fit_single(
+        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting):
     """
     Fit a ramp using ordinary least squares. Calculate the count rate for each
     pixel in all data cube sections and all integrations, equal to the weighted
@@ -645,6 +655,9 @@ def ols_ramp_fit_single(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, we
     weighting : str
         'optimal' is the only valid value
 
+    avg_dark_current : float (electrons)
+        The average dark current for this detector
+
     Return
     ------
     image_info : tuple
@@ -1049,6 +1062,8 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans)
     fit_slopes_ans : tuple
         Contains intermediate values computed in the first pass over the data.
 
+    avg_dark_current : float (electrons)
+        The average dark current for this detector
     Return
     ------
     var_p3 : ndarray
@@ -1138,7 +1153,8 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans)
             # Suppress harmless arithmetic warnings for now
             warnings.filterwarnings("ignore", ".*invalid value.*", RuntimeWarning)
             warnings.filterwarnings("ignore", ".*divide by zero.*", RuntimeWarning)
-            var_p4[num_int, :, rlo:rhi, :] = den_p3 * med_rates[rlo:rhi, :]
+            var_p4[num_int, :, rlo:rhi, :] = ((den_p3 * med_rates[rlo:rhi, :]) +
+                                              ramp_data.avg_dark_current)
 
             # Find the segment variance due to read noise and convert back to DN
             var_r4[num_int, :, rlo:rhi, :] = num_r3 * den_r3 / gain_sect**2
@@ -1178,7 +1194,10 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans)
         var_p3[:, med_rates <= 0.0] = 0.0
         warnings.resetwarnings()
 
-        var_p4[num_int, :, med_rates <= 0.0] = 0.0
+        # For pixels with zero or negative rates set the Poisson variance to the
+        # dark current rate.
+        var_p4[num_int, :, med_rates <= 0.0] = ramp_data.avg_dark_current
+
         var_both4[num_int, :, :, :] = var_r4[num_int, :, :, :] + var_p4[num_int, :, :, :]
         inv_var_both4[num_int, :, :, :] = 1.0 / var_both4[num_int, :, :, :]
 
diff --git a/src/stcal/ramp_fitting/ramp_fit.py b/src/stcal/ramp_fitting/ramp_fit.py
index 514429f1b..554c77674 100755
--- a/src/stcal/ramp_fitting/ramp_fit.py
+++ b/src/stcal/ramp_fitting/ramp_fit.py
@@ -30,7 +30,7 @@
 BUFSIZE = 1024 * 300000  # 300Mb cache size for data section
 
 
-def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False):
+def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False, avg_dark_current=0):
     """
     Create an internal ramp fit class from a data model.
 
@@ -77,6 +77,7 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False):
     ramp_data.num_rows = ramp_data.data.shape[2]
 
     ramp_data.suppress_one_group_ramps = suppress_one_group
+    ramp_data.avg_dark_current = avg_dark_current
 
     return ramp_data
 
@@ -92,7 +93,8 @@ def ramp_fit(
     max_cores,
     dqflags,
     suppress_one_group=False,
-):
+    avg_dark_current=0.0,
+    ):
     """
     Calculate the count rate for each pixel in all data cube sections and all
     integrations, equal to the slope for all sections (intervals between
@@ -141,6 +143,9 @@ def ramp_fit(
         Find ramps with only one good group and treat it like it has zero good
         groups.
 
+    avg_dark_current : flt
+        The average dark current for this detector in units of electrons per second.
+
     Returns
     -------
     image_info : tuple
@@ -164,16 +169,15 @@ def ramp_fit(
     # Create an instance of the internal ramp class, using only values needed
     # for ramp fitting from the to remove further ramp fitting dependence on
     # data models.
-    ramp_data = create_ramp_fit_class(model, dqflags, suppress_one_group)
+    ramp_data = create_ramp_fit_class(model, dqflags, suppress_one_group, avg_dark_current)
 
     return ramp_fit_data(
-        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, algorithm, weighting, max_cores, dqflags
-    )
+        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d,
+        algorithm, weighting, max_cores, dqflags, avg_dark_current)
 
 
-def ramp_fit_data(
-    ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, algorithm, weighting, max_cores, dqflags
-):
+def ramp_fit_data(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d,
+                  algorithm, weighting, max_cores):
     """
     This function begins the ramp fit computation after the creation of the
     RampData class.  It determines the proper path for computation to take
@@ -211,10 +215,6 @@ def ramp_fit_data(
         to use for multi-proc. The total number of cores includes the SMT cores
         (Hyper Threading for Intel).
 
-    dqflags : dict
-        A dictionary with at least the following keywords:
-        DO_NOT_USE, SATURATED, JUMP_DET, NO_GAIN_VALUE, UNRELIABLE_SLOPE
-
     Returns
     -------
     image_info : tuple
@@ -247,8 +247,7 @@ def ramp_fit_data(
 
         # Compute ramp fitting using ordinary least squares.
         image_info, integ_info, opt_info = ols_fit.ols_ramp_fit_multi(
-            ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, max_cores
-        )
+            ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, max_cores)
         gls_opt_info = None
 
     return image_info, integ_info, opt_info, gls_opt_info
diff --git a/src/stcal/ramp_fitting/ramp_fit_class.py b/src/stcal/ramp_fitting/ramp_fit_class.py
index f8a78efd2..8fd089f89 100644
--- a/src/stcal/ramp_fitting/ramp_fit_class.py
+++ b/src/stcal/ramp_fitting/ramp_fit_class.py
@@ -15,6 +15,7 @@ def __init__(self):
         self.groupgap = None
         self.nframes = None
         self.drop_frames1 = None
+        self.avg_dark_current = 0.0
 
         # Data quality flags
         self.flags_do_not_use = None
diff --git a/tests/test_jump.py b/tests/test_jump.py
index a85eef519..36bea8198 100644
--- a/tests/test_jump.py
+++ b/tests/test_jump.py
@@ -1,6 +1,6 @@
 import numpy as np
 import pytest
-
+from astropy.io import fits
 from stcal.jump.jump import (
     calc_num_slices,
     extend_saturation,
@@ -11,6 +11,7 @@
     detect_jumps,
 )
 
+
 DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8}
 
 
@@ -445,4 +446,26 @@ def test_calc_num_slices():
     assert calc_num_slices(n_rows, "all", max_available_cores) == 10
     assert calc_num_slices(n_rows, "3/4", max_available_cores) == 1
     n_rows = 9
-    assert calc_num_slices(n_rows, "21", max_available_cores) == 9
+    assert (calc_num_slices(n_rows, '21', max_available_cores) == 9)
+
+@pytest.mark.skip(reason="local test only")
+def test_shower_slowmode():
+    data = fits.getdata("shower_input_data.fits")
+    gdq = fits.getdata("shower_input_gdq.fits")
+    nint = data.shape[0]
+    ngrps = data.shape[1]
+    ncols = data.shape[2]
+    nrows = data.shape[3]
+    gain = 5
+    readnoise = fits.getdata("/grp/crds/jwst/references/jwst/jwst_miri_readnoise_0086.fits")
+#    readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain
+    rng = np.random.default_rng(12345)
+#    data[0, 1:, 14:20, 15:20] = 6 * gain * 1.7
+#    data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise
+    gdq, num_showers = find_faint_extended(data, gdq, readnoise, 1, 100,
+                                           snr_threshold=1.3,
+                                           min_shower_area=40, inner=1,
+                                           outer=2, sat_flag=2, jump_flag=4,
+                                           ellipse_expand=1.1, num_grps_masked=3)
+    print("number of showers", num_showers)
+    fits.writeto("outgdq.fits", gdq, overwrite=True)
\ No newline at end of file
diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py
index d8e906104..ccb35c7f4 100644
--- a/tests/test_ramp_fitting.py
+++ b/tests/test_ramp_fitting.py
@@ -3,6 +3,7 @@
 from stcal.ramp_fitting.ramp_fit import ramp_fit_data
 from stcal.ramp_fitting.ramp_fit_class import RampData
 from stcal.ramp_fitting.utils import compute_num_slices
+from stcal.ramp_fitting.ols_fit import ramp_fit_compute_variances
 
 DELIM = "=" * 70
 
@@ -41,7 +42,7 @@ def base_neg_med_rates_single_integration():
     var = (rnoise_val, gain_val)
     tm = (nframes, gtime, dtime)
     ramp_data, rnoise, gain = setup_inputs(dims, var, tm)
-
+    avg_dark_current = 0.0
     # Set up negative ramp
     neg_ramp = np.array([k + 1 for k in range(ngroups)])
     nslope = -0.5
@@ -51,7 +52,7 @@ def base_neg_med_rates_single_integration():
     # Run ramp fit on RampData
     buffsize, save_opt, algo, wt, ncores = 512, True, "OLS", "optimal", "none"
     slopes, cube, optional, gls_dummy = ramp_fit_data(
-        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores, dqflags
+        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores,
     )
 
     return slopes, cube, optional, gls_dummy
@@ -113,6 +114,7 @@ def base_neg_med_rates_multi_integrations():
     nints, ngroups, nrows, ncols = 3, 10, 1, 1
     rnoise_val, gain_val = 10.0, 1.0
     nframes, gtime, dtime = 1, 1.0, 1
+    avg_dark_current = 0.0
     dims = (nints, ngroups, nrows, ncols)
     var = (rnoise_val, gain_val)
     tm = (nframes, gtime, dtime)
@@ -130,7 +132,7 @@ def base_neg_med_rates_multi_integrations():
     # Run ramp fit on RampData
     buffsize, save_opt, algo, wt, ncores = 512, True, "OLS", "optimal", "none"
     slopes, cube, optional, gls_dummy = ramp_fit_data(
-        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores, dqflags
+        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores
     )
 
     return slopes, cube, optional, gls_dummy, dims
@@ -198,6 +200,7 @@ def base_neg_med_rates_single_integration_multi_segment():
     var = (rnoise_val, gain_val)
     tm = (nframes, gtime, dtime)
     ramp_data, rnoise, gain = setup_inputs(dims, var, tm)
+    avg_dark_current = 0.0
 
     # Set up negative ramp
     neg_ramp = np.array([k + 1 for k in range(ngroups)])
@@ -213,7 +216,7 @@ def base_neg_med_rates_single_integration_multi_segment():
     # Run ramp fit on RampData
     buffsize, save_opt, algo, wt, ncores = 512, True, "OLS", "optimal", "none"
     slopes, cube, optional, gls_dummy = ramp_fit_data(
-        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores, dqflags
+        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores,
     )
 
     return slopes, cube, optional, gls_dummy, dims
@@ -260,6 +263,7 @@ def test_utils_dq_compress_final():
     nints, ngroups, nrows, ncols = 2, 5, 1, 3
     rnoise_val, gain_val = 10.0, 1.0
     nframes, gtime, dtime = 1, 1.0, 1
+    avg_dark_current = 0.0
     dims = (nints, ngroups, nrows, ncols)
     var = (rnoise_val, gain_val)
     tm = (nframes, gtime, dtime)
@@ -273,7 +277,7 @@ def test_utils_dq_compress_final():
     # Run ramp fit on RampData
     buffsize, save_opt, algo, wt, ncores = 512, True, "OLS", "optimal", "none"
     slopes, cube, optional, gls_dummy = ramp_fit_data(
-        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores, dqflags
+        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores,
     )
 
     dq = slopes[1]
@@ -316,7 +320,7 @@ def jp_2326_test_setup():
     )
     dnu = dqflags["DO_NOT_USE"]
     dq = np.array([dnu, 0, 0, 0, 0, 0, 0, 0, 0, dnu])
-
+    avg_dark_current = 0.0
     nints, ngroups, nrows, ncols = 1, len(ramp), 1, 1
     data = np.zeros((nints, ngroups, nrows, ncols))
     gdq = np.zeros((nints, ngroups, nrows, ncols), dtype=np.uint8)
@@ -329,9 +333,10 @@ def jp_2326_test_setup():
     ramp_data = RampData()
     ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pdq)
     ramp_data.set_meta(
-        name="MIRI", frame_time=2.77504, group_time=2.77504, groupgap=0, nframes=1, drop_frames1=None
-    )
+        name="MIRI", frame_time=2.77504, group_time=2.77504, groupgap=0,
+        nframes=1, drop_frames1=None)
     ramp_data.set_dqflags(dqflags)
+    ramp_data.avg_dark_current = avg_dark_current
 
     # Set up gain and read noise
     gain = np.ones(shape=(nrows, ncols), dtype=np.float32) * 5.5
@@ -347,11 +352,11 @@ def test_miri_ramp_dnu_at_ramp_beginning():
     """
     ramp_data, gain, rnoise = jp_2326_test_setup()
     ramp_data.groupdq[0, 1, 0, 0] = dqflags["DO_NOT_USE"]
-
+    avg_dark_current = 0.0
     # Run ramp fit on RampData
     buffsize, save_opt, algo, wt, ncores = 512, True, "OLS", "optimal", "none"
     slopes1, cube, optional, gls_dummy = ramp_fit_data(
-        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores, dqflags
+        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores,
     )
 
     s1 = slopes1[0]
@@ -369,11 +374,11 @@ def test_miri_ramp_dnu_and_jump_at_ramp_beginning():
     """
     ramp_data, gain, rnoise = jp_2326_test_setup()
     ramp_data.groupdq[0, 1, 0, 0] = dqflags["JUMP_DET"]
-
+    avg_dark_current = 0.0
     # Run ramp fit on RampData
     buffsize, save_opt, algo, wt, ncores = 512, True, "OLS", "optimal", "none"
     slopes2, cube, optional, gls_dummy = ramp_fit_data(
-        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores, dqflags
+        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores,
     )
 
     s2 = slopes2[0]
@@ -392,7 +397,7 @@ def test_2_group_cases():
     base_err = [0.0, 0.0]
     gain_val = 0.9699
     rnoise_val = 9.4552
-
+    avg_dark_current = 0.0
     possibilities = [
         # Both groups are good
         [GOOD, GOOD],
@@ -444,7 +449,7 @@ def test_2_group_cases():
     # Run ramp fit on RampData
     buffsize, save_opt, algo, wt, ncores = 512, True, "OLS", "optimal", "none"
     slopes, cube, optional, gls_dummy = ramp_fit_data(
-        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores, dqflags
+        ramp_data, buffsize, save_opt, rnoise, gain, algo, wt, ncores,
     )
 
     # Check the outputs
@@ -515,11 +520,12 @@ def run_one_group_ramp_suppression(nints, suppress):
         ramp_data.groupdq[1, :, 0, 2] = zdq
 
     ramp_data.suppress_one_group_ramps = suppress
-
+    avg_dark_current = 0.0
     algo = "OLS"
     save_opt, ncores, bufsize = False, "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo,
+        "optimal", ncores,
     )
 
     return slopes, cube, dims
@@ -753,12 +759,11 @@ def create_zero_frame_data():
     # Set up group DQ array.
     gdq[0, :, :, :] = dqflags["SATURATED"]
     gdq[0, 0, 0, 0] = dqflags["GOOD"]
-
+    avg_dark_current = 0.0
     # Create RampData for testing.
     ramp_data = RampData()
     ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq)
-    ramp_data.set_meta(
-        name="NIRCam",
+    ramp_data.set_meta(name="NIRCam",
         frame_time=frame_time,
         group_time=group_time,
         groupgap=groupgap,
@@ -767,6 +772,7 @@ def create_zero_frame_data():
     )
     ramp_data.set_dqflags(dqflags)
 
+
     ramp_data.suppress_one_group_ramps = False
     ramp_data.zeroframe = zframe
 
@@ -789,11 +795,11 @@ def test_zeroframe():
     The second integration has all good groups with half the data values.
     """
     ramp_data, gain, rnoise = create_zero_frame_data()
-
+    avg_dark_current = 0.0
     algo, save_opt, ncores, bufsize = "OLS", False, "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
-    )
+        ramp_data, bufsize, save_opt, rnoise, gain, algo,
+        "optimal", ncores,)
 
     tol = 1.0e-5
 
@@ -874,7 +880,7 @@ def create_only_good_0th_group_data():
 
     gdq[0, :, 0, 2] = np.array([SAT] * ngroups)
     gdq[0, 0, 0, 2] = GOOD
-
+    avg_dark_current = 0.0
     # Create RampData for testing.
     ramp_data = RampData()
     ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq)
@@ -907,10 +913,11 @@ def test_only_good_0th_group():
 
     # Dimensions are (1, 5, 1, 3)
     ramp_data, gain, rnoise = create_only_good_0th_group_data()
-
+    avg_dark_current = 0.0
     algo, save_opt, ncores, bufsize = "OLS", False, "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo,
+        "optimal", ncores,
     )
 
     tol = 1.0e-5
@@ -956,10 +963,11 @@ def test_all_sat():
 
     ramp, gain, rnoise = create_blank_ramp_data(dims, var, tm)
     ramp.groupdq[:, 0, :, :] = ramp.flags_saturated
-
+    avg_dark_current = 0.0
     algo, save_opt, ncores, bufsize = "OLS", False, "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp, bufsize, save_opt, rnoise, gain, algo,
+        "optimal", ncores,
     )
 
     assert slopes is None
@@ -986,11 +994,10 @@ def test_dq_multi_int_dnu():
     ramp.data[0, :, 0, 0] = np.array(base_arr)
     ramp.data[1, :, 0, 0] = np.array(base_arr)
     ramp.groupdq[0, :, 0, 0] = np.array(dq_arr)
-
+    avg_dark_current = 0.0
     algo, save_opt, ncores, bufsize = "OLS", False, "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
-    )
+        ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,)
 
     tol = 1.0e-5
 
@@ -1078,6 +1085,7 @@ def test_multi_more_cores_than_rows():
         ]
     )
     factor = 1.05
+    avg_dark_current = 0.0
     for integ in range(nints):
         for row in range(nrows):
             for col in range(ncols):
@@ -1086,8 +1094,8 @@ def test_multi_more_cores_than_rows():
 
     bufsize, algo, save_opt, ncores = 512, "OLS", False, "all"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
-    )
+        ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
+         )
     # This part of the test is simply to make sure ramp fitting
     # doesn't crash.  No asserts are necessary here.
 
@@ -1234,10 +1242,10 @@ def test_new_saturation():
     use.
     """
     ramp, gain, rnoise = get_new_saturation()
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = False, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp, bufsize, save_opt, rnoise, gain, algo,"optimal", ncores,
     )
 
     tol = 1.0e-5
@@ -1322,12 +1330,12 @@ def test_invalid_integrations():
     for integ in range(nints):
         ramp.data[integ, :, 0, 0] = np.array(int_data[integ], dtype=np.float32)
         ramp.groupdq[integ, :, 0, 0] = np.array(int_dq[integ], dtype=np.uint8)
-
+    avg_dark_current = 0.0
     ramp.suppress_one_group_ramps = True
 
     save_opt, ncores, bufsize, algo = False, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp, bufsize, save_opt, rnoise, gain, algo,"optimal", ncores,
     )
 
     tol = 1.0e-5
@@ -1386,10 +1394,10 @@ def test_one_group():
     ramp, gain, rnoise = create_blank_ramp_data(dims, var, tm)
 
     ramp.data[0, 0, 0, 0] = 105.31459
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = False, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp, bufsize, save_opt, rnoise, gain, algo,"optimal", ncores,
     )
 
     tol = 1e-5
@@ -1458,6 +1466,28 @@ def test_compute_num_slices():
     n_rows = 9
     assert compute_num_slices("21", n_rows, max_available_cores) == 9
 
+def test_ramp_fit_dark_current_variance():
+    nints, ngroups, nrows, ncols = 1, 5, 2, 2
+    rnval, gval = 10.0, 5.0
+    frame_time, nframes, groupgap = 10, 1, 0
+    var = rnval, gval
+    tm = frame_time, nframes, groupgap
+    dims = nints, ngroups, nrows, ncols
+    ramp, gain, rnoise = create_blank_ramp_data(dims, var, tm)
+    #Test for zero dark rate
+    ramp.avg_dark_current = 0.0
+    med_rates = np.zeros(shape=(2, 2))
+    fit_slopes_ans = (1, 0, 0, 0, 0, 1, 0, 0, 0, med_rates)
+    var_p3, var_r3, var_p4, var_r4, var_both4, var_both3, inv_var_both4, s_inv_var_p3, \
+    s_inv_var_r3, s_inv_var_both3 = ramp_fit_compute_variances(ramp, gain, rnoise, fit_slopes_ans)
+
+    assert(var_p4[0, 0, 0, 0] == 0.0)
+
+    #Test for a dark rate of 1.5 e-/sec
+    ramp.avg_dark_current = 1.5
+    var_p3, var_r3, var_p4, var_r4, var_both4, var_both3, inv_var_both4, s_inv_var_p3, \
+        s_inv_var_r3, s_inv_var_both3 = ramp_fit_compute_variances(ramp, gain, rnoise, fit_slopes_ans)
+    assert(var_p4[0, 0, 0, 0] == 1.5)
 
 # -----------------------------------------------------------------------------
 #                           Set up functions
@@ -1486,12 +1516,12 @@ def setup_inputs(dims, var, tm):
 
     for c_int in range(1, nints):
         data[c_int, :, :, :] = data[0, :, :, :].copy()
-
+    avg_dark_current = 0.0
     ramp_data = RampData()
     ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq)
     ramp_data.set_meta(
-        name="MIRI", frame_time=dtime, group_time=gtime, groupgap=0, nframes=nframes, drop_frames1=None
-    )
+        name="MIRI", frame_time=dtime, group_time=gtime, groupgap=0,
+        nframes=nframes, drop_frames1=None)
     ramp_data.set_dqflags(dqflags)
 
     gain = np.ones(shape=(nrows, ncols), dtype=np.float64) * gain
diff --git a/tests/test_ramp_fitting_cases.py b/tests/test_ramp_fitting_cases.py
index 675e6a759..11fef3f8e 100644
--- a/tests/test_ramp_fitting_cases.py
+++ b/tests/test_ramp_fitting_cases.py
@@ -57,10 +57,10 @@ def test_pix_0():
     )
     dq = [GOOD] * ngroups
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -98,10 +98,10 @@ def test_pix_1():
     dq[2] = JUMP
     dq[4:] = [SAT] * (ngroups - 4)
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -131,10 +131,10 @@ def test_pix_2():
     )
     dq = [GOOD, GOOD, GOOD, JUMP, GOOD, JUMP, GOOD, JUMP, SAT, SAT]
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -174,10 +174,10 @@ def test_pix_3():
     dq = [GOOD] * ngroups
     dq[-2] = JUMP
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -216,10 +216,10 @@ def test_pix_4():
     )
     dq = [GOOD] + [SAT] * (ngroups - 1)
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -299,10 +299,10 @@ def test_pix_5():
     dq = [GOOD] * ngroups
     dq[4] = JUMP
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # XXX see the note above for the differences in C and python testing values.
@@ -348,10 +348,10 @@ def test_pix_6():
     dq[2] = JUMP
     dq[3] = JUMP
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -389,10 +389,10 @@ def test_pix_7():
     )
     dq = [GOOD] * (ngroups - 2) + [JUMP, JUMP]
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -422,10 +422,10 @@ def test_pix_8():
     )
     dq = [GOOD, JUMP, GOOD, GOOD, GOOD, GOOD, GOOD, SAT, SAT, SAT]
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -456,10 +456,10 @@ def test_pix_9():
     )
     dq = [GOOD, GOOD, JUMP, JUMP, GOOD, GOOD, GOOD, GOOD, JUMP, GOOD]
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -499,10 +499,10 @@ def test_pix_10():
     )
     dq = [GOOD, GOOD, JUMP, GOOD, GOOD, JUMP, GOOD, GOOD, GOOD, GOOD]
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -540,10 +540,10 @@ def test_pix_11():
     )
     dq = [GOOD, GOOD] + [SAT] * (ngroups - 2)
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -576,10 +576,10 @@ def test_pix_12():
     ramp_data.groupdq[0, :, 0, 0] = np.array([0, SAT])
     ramp_data.data[0, :, 0, 1] = np.array([61000.0, 61000.0], dtype=np.float32)
     ramp_data.groupdq[0, :, 0, 1] = np.array([SAT, SAT])
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results for pixel 1:
@@ -630,10 +630,10 @@ def test_miri_0():
     )
     dq = [DNU] + [GOOD] * (ngroups - 2) + [DNU]
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -664,10 +664,10 @@ def test_miri_1():
     )
     dq = [DNU | JUMP] + [GOOD] * (ngroups - 2) + [DNU]
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -698,10 +698,10 @@ def test_miri_2():
     )
     dq = [DNU | JUMP] + [GOOD] * (ngroups - 2) + [DNU | JUMP]
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
@@ -732,10 +732,10 @@ def test_miri_3():
     )
     dq = [DNU] + [GOOD] * (ngroups - 2) + [DNU | JUMP]
     ramp_data.groupdq[0, :, 0, 0] = np.array(dq)
-
+    avg_dark_current = 0.0
     save_opt, ncores, bufsize, algo = True, "none", 1024 * 30000, "OLS"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores, dqflags
+        ramp_data, bufsize, save_opt, rnoise, gain, algo, "optimal", ncores,
     )
 
     # Set truth values for PRIMARY results:
diff --git a/tests/test_ramp_fitting_gls_fit.py b/tests/test_ramp_fitting_gls_fit.py
index e940bd949..2076d957a 100644
--- a/tests/test_ramp_fitting_gls_fit.py
+++ b/tests/test_ramp_fitting_gls_fit.py
@@ -103,10 +103,10 @@ def test_one_group_small_buffer():
     ramp_data, gain2d, rnoise2d = setup_inputs(dims, gain, rnoise, group_time, frame_time)
 
     ramp_data.data[0, 0, 50, 50] = 10.0
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores = False, "GLS", "none"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores,
     )
 
     data = slopes[0]
@@ -130,10 +130,10 @@ def test_two_integrations():
     ramp = np.asarray([x * 100 for x in range(11)])
     ramp_data.data[0, :, row, col] = ramp
     ramp_data.data[1, :, row, col] = ramp * 2
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores = False, "GLS", "none"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores,
     )
 
     np.testing.assert_allclose(slopes[0][row, col], 133.3377685, 1e-6)
@@ -152,10 +152,10 @@ def test_one_group_two_integrations():
 
     ramp_data.data[0, 0, 0, 0] = 10.0
     ramp_data.data[1, 0, 0, 0] = 11.0
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores = False, "GLS", "none"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     data = slopes[0]
@@ -174,10 +174,10 @@ def test_nocrs_noflux():
     group_time, frame_time = 1.0, 1
 
     ramp_data, gain2d, rnoise2d = setup_inputs(dims, gain, rnoise, group_time, frame_time)
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores = False, "GLS", "none"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     assert np.max(slopes[0]) == 0
@@ -197,10 +197,10 @@ def test_nocrs_noflux_firstrows_are_nan():
 
     ramp_data, gain2d, rnoise2d = setup_inputs(dims, gain, rnoise, group_time, frame_time)
     ramp_data.data[0:, 0:12, :] = np.nan
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores = False, "GLS", "none"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     assert np.max(slopes[0]) == 0
@@ -220,11 +220,11 @@ def test_error_when_frame_time_not_set():
     ramp_data.data[0:, 0:12, :] = np.nan
 
     ramp_data.frame_time = None  # Must be set
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores = False, "GLS", "none"
     with pytest.raises(UnboundLocalError):
         slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-            ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+            ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
         )
 
 
@@ -250,10 +250,10 @@ def test_five_groups_two_integrations_Poisson_noise_only():
     ramp_data.data[1, 2, row, col] = 25.0
     ramp_data.data[1, 3, row, col] = 33.0
     ramp_data.data[1, 4, row, col] = 160.0
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores = False, "GLS", "none"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     out_slope = slopes[0][row, col]
@@ -278,11 +278,11 @@ def test_bad_gain_values():
     ramp_data, gain2d, rnoise2d = setup_inputs(dims, gain, rnoise, group_time, frame_time)
     gain2d[r1, c1] = -10
     gain2d[r2, c2] = np.nan
-
+    avg_dark_current = 0.0
     # save_opt, algo, ncores = False, "OLS", "none"
     save_opt, algo, ncores = False, "GLS", "none"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     # data, dq, var_poisson, var_rnoise, err = slopes
@@ -311,10 +311,10 @@ def test_simple_ramp():
 
     ramp = np.array(list(range(ngroups))) * 20 + 10
     ramp_data.data[0, :, 50, 50] = ramp
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores = False, "GLS", "none"
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, 512, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -336,10 +336,10 @@ def test_read_noise_only_fit():
 
     ramp_arr = [10.0, 15.0, 25.0, 33.0, 60.0]
     ramp_data.data[0, :, 50, 50] = np.array(ramp_arr)
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     xvalues = np.arange(5) * 1.0
@@ -367,10 +367,10 @@ def test_photon_noise_only_fit():
     ramp_data.data[0, :, 50, 50] = np.array(ramp_arr)
 
     check = (ramp_data.data[0, 4, 50, 50] - ramp_data.data[0, 0, 50, 50]) / 4.0
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -398,10 +398,10 @@ def test_photon_noise_only_bad_last_group():
     check = (ramp_data.data[0, 3, 50, 50] - ramp_data.data[0, 0, 50, 50]) / 3.0
 
     ramp_data.groupdq[0, 4, :, :] = DO_NOT_USE
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -426,10 +426,10 @@ def test_photon_noise_with_unweighted_fit():
     ramp_data.data[0, 2, 50, 50] = 25.0
     ramp_data.data[0, 3, 50, 50] = 33.0
     ramp_data.data[0, 4, 50, 50] = 60.0
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "unweighted", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "unweighted", ncores, 
     )
 
     xvalues = np.arange(5) * 1.0
@@ -468,10 +468,10 @@ def test_two_groups_fit():
     # 1st group is saturated
     ramp_data.groupdq[0, 0, 0, 2] = SATURATED
     ramp_data.groupdq[0, 1, 0, 2] = SATURATED
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     ans_data = slopes[0][0, 0]
@@ -500,10 +500,10 @@ def test_four_groups_oneCR_orphangroupatend_fit():
     ramp_data.groupdq[0, 3, 50, 50] = JUMP_DET
 
     check = ramp_data.data[0, 1, 50, 50] - ramp_data.data[0, 0, 50, 50]
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -529,10 +529,10 @@ def test_four_groups_two_CRs_at_end():
 
     ramp_data.groupdq[0, 2, 50, 50] = JUMP_DET
     ramp_data.groupdq[0, 3, 50, 50] = JUMP_DET
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -558,10 +558,10 @@ def test_four_groups_four_CRs():
     ramp_data.groupdq[0, 1, 50, 50] = JUMP_DET
     ramp_data.groupdq[0, 2, 50, 50] = JUMP_DET
     ramp_data.groupdq[0, 3, 50, 50] = JUMP_DET
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -590,10 +590,10 @@ def test_four_groups_three_CRs_at_end():
     ramp_data.groupdq[0, 1, 50, 50] = JUMP_DET
     ramp_data.groupdq[0, 2, 50, 50] = JUMP_DET
     ramp_data.groupdq[0, 3, 50, 50] = JUMP_DET
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -619,10 +619,10 @@ def test_four_groups_CR_causes_orphan_1st_group():
     ramp_data.data[0, 3, 50, 50] = 165.0
 
     ramp_data.groupdq[0, 1, 50, 50] = JUMP_DET
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -640,10 +640,10 @@ def test_one_group_fit():
 
     ramp_data, gain2d, rnoise2d = setup_inputs(dims, gain, rnoise, group_time, frame_time)
     ramp_data.data[0, 0, 50, 50] = 10.0
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -665,10 +665,10 @@ def test_two_groups_unc():
 
     ramp_data.data[0, 0, 50, 50] = 10.0
     ramp_data.data[0, 1, 50, 50] = 10.0 + deltaDN
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[2][50, 50]
@@ -705,10 +705,10 @@ def test_five_groups_unc():
     ramp_data.data[0, 4, 50, 50] = 60.0
 
     # check = np.median(np.diff(ramp_data.data[0,:,50,50])) / grouptime
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     """
@@ -752,10 +752,10 @@ def test_oneCR_10_groups_combination():
     ramp_data.data[0, 9, 50, 50] = 180.0
 
     ramp_data.groupdq[0, 5, 50, 50] = JUMP_DET
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, 
     )
 
     answer = slopes[0][50, 50]
@@ -806,10 +806,10 @@ def test_oneCR_10_groups_combination_noisy2ndSegment():
     ramp_data.data[0, 9, 50, 50] = 180.0
 
     ramp_data.groupdq[0, 5, 50, 50] = JUMP_DET
-
+    avg_dark_current = 0.0
     save_opt, algo, ncores, bufsize = False, "GLS", "none", 1024 * 30000
     slopes, cube, ols_opt, gls_opt = ramp_fit_data(
-        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores, test_dq_flags
+        ramp_data, bufsize, save_opt, rnoise2d, gain2d, algo, "optimal", ncores,
     )
 
     """