diff --git a/CHANGES.rst b/CHANGES.rst
index c1eb217f4..a01b689ba 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -7,7 +7,8 @@ general
 
 - DQ step flags science data affected by guide window read [#599]
 
-
+- Added support for Quantities for data arrays. [#613]
+  
 0.9.0 (2022-11-14)
 ==================
 
diff --git a/pyproject.toml b/pyproject.toml
index 37494f222..42ea23b3a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,7 +12,7 @@ classifiers = [
     'Programming Language :: Python :: 3',
 ]
 dependencies = [
-    'asdf >=2.12.1',
+    'asdf >=2.13.0',
     'astropy >=5.0.4',
     'crds >=11.16.16',
     'gwcs >=0.18.1',
@@ -20,8 +20,9 @@ dependencies = [
     'numpy >=1.20',
     'pyparsing >=2.4.7',
     'requests >=2.22',
-    'roman_datamodels >=0.14.0',
+    #'roman_datamodels >=0.14.0',
     #'roman_datamodels @ git+https://github.com/spacetelescope/roman_datamodels.git@main',
+    'roman_datamodels @ git+https://github.com/PaulHuwe/roman_datamodels.git@ArrayQuanitites',
     'stcal >=1.2.1, <2.0',
     'stpipe >=0.4.2, <1.0',
 ]
diff --git a/romancal/dark_current/dark_current_step.py b/romancal/dark_current/dark_current_step.py
index b51d9a1a6..cd7e27202 100755
--- a/romancal/dark_current/dark_current_step.py
+++ b/romancal/dark_current/dark_current_step.py
@@ -5,6 +5,8 @@
 from roman_datamodels import datamodels as rdd
 from stcal.dark_current import dark_sub
 from roman_datamodels.testing import utils as testutil
+from astropy import units as u
+from roman_datamodels import units as ru
 
 
 __all__ = ["DarkCurrentStep"]
@@ -43,7 +45,7 @@ def process(self, input):
                 dark_model.meta.exposure['groupgap'] = input_model.meta.exposure.groupgap
 
             # Reshaping data variables for stcal compatibility
-            input_model.data = input_model.data.astype(np.float32)[np.newaxis, :]
+            input_model.data = input_model.data[np.newaxis, :]
             input_model.groupdq = input_model.groupdq[np.newaxis, :]
             input_model.err = input_model.err[np.newaxis, :]
 
@@ -133,10 +135,10 @@ def dark_output_data_as_ramp_model(out_data, input_model):
     # Removing integration dimension from variables (added for stcal
     # compatibility)
     # Roman 3D
-    out_model.data = out_data.data[0]
+    out_model.data = u.Quantity(out_data.data[0], ru.DN, dtype=out_data.data.dtype)
     out_model.groupdq = out_data.groupdq[0]
     # Roman 2D
     out_model.pixeldq = out_data.pixeldq
-    out_model.err = out_data.err[0]
+    out_model.err = u.Quantity(out_data.err[0], ru.DN, dtype=out_data.err.dtype)
 
     return out_model
diff --git a/romancal/dark_current/tests/test_dark.py b/romancal/dark_current/tests/test_dark.py
index 709491b47..acfdc43dc 100755
--- a/romancal/dark_current/tests/test_dark.py
+++ b/romancal/dark_current/tests/test_dark.py
@@ -5,12 +5,15 @@
 import pytest
 import numpy as np
 import os
+from astropy import units as u
 
 from romancal.dark_current import DarkCurrentStep
 
 from roman_datamodels.datamodels import RampModel, DarkRefModel
 import roman_datamodels as rdm
 from roman_datamodels.testing import utils as testutil
+from roman_datamodels import units as ru
+
 
 
 @pytest.mark.parametrize(
@@ -71,17 +74,17 @@ def test_dark_step_subtraction(instrument, exptype):
 
     # populate data array of science cube
     for i in range(0, 20):
-        ramp_model.data[0, 0, i] = i
+        ramp_model.data.value[0, 0, i] = i
         darkref_model.data[0, 0, i] = i * 0.1
 
     # Perform Dark Current subtraction step
     result = DarkCurrentStep.call(ramp_model, override_dark=darkref_model)
 
     # check that the dark file is subtracted frame by frame from the science data
-    diff = ramp_model.data - darkref_model.data
+    diff = ramp_model.data.value - darkref_model.data
 
     # test that the output data file is equal to the difference found when subtracting reffile from sci file
-    np.testing.assert_array_equal(result.data, diff, err_msg='dark file should be subtracted from sci file ')
+    np.testing.assert_array_equal(result.data.value, diff, err_msg='dark file should be subtracted from sci file ')
 
 
 @pytest.mark.parametrize(
@@ -127,7 +130,7 @@ def create_ramp_and_dark(shape, instrument, exptype):
     ramp.meta.instrument.detector = 'WFI01'
     ramp.meta.instrument.optical_element = 'F158'
     ramp.meta.exposure.type = exptype
-    ramp.data = np.ones(shape, dtype=np.float32)
+    ramp.data = u.Quantity(np.ones(shape, dtype=np.float32) , ru.DN, dtype=np.float32)
     ramp_model = RampModel(ramp)
 
     # Create dark model
diff --git a/romancal/dq_init/tests/test_dq_init.py b/romancal/dq_init/tests/test_dq_init.py
index 9d4381a49..8b3eccfef 100644
--- a/romancal/dq_init/tests/test_dq_init.py
+++ b/romancal/dq_init/tests/test_dq_init.py
@@ -3,10 +3,12 @@
 import numpy as np
 import pytest
 import warnings
+from astropy import units as u
 
 from roman_datamodels import stnode
 from roman_datamodels.datamodels import MaskRefModel, ScienceRawModel
 from roman_datamodels.testing import utils as testutil
+from roman_datamodels import units as ru
 
 from romancal.lib import dqflags
 from romancal.dq_init import DQInitStep
@@ -196,7 +198,7 @@ def test_dqinit_step_interface(instrument, exptype):
     wfi_sci_raw.meta['guidestar']['gw_window_xstart'] = 1012
     wfi_sci_raw.meta['guidestar']['gw_window_xsize'] = 16
     wfi_sci_raw.meta.exposure.type = exptype
-    wfi_sci_raw.data = np.ones(shape, dtype=np.uint16)
+    wfi_sci_raw.data = u.Quantity(np.ones(shape, dtype=np.uint16), ru.DN, dtype=np.uint16)
     wfi_sci_raw_model = ScienceRawModel(wfi_sci_raw)
 
     # Create mask model
@@ -248,7 +250,7 @@ def test_dqinit_refpix(instrument, exptype):
     wfi_sci_raw.meta['guidestar']['gw_window_xstart'] = 1012
     wfi_sci_raw.meta['guidestar']['gw_window_xsize'] = 16
     wfi_sci_raw.meta.exposure.type = exptype
-    wfi_sci_raw.data = np.ones(shape, dtype=np.uint16)
+    wfi_sci_raw.data = u.Quantity(np.ones(shape, dtype=np.uint16), ru.DN, dtype=np.uint16)
     wfi_sci_raw_model = ScienceRawModel(wfi_sci_raw)
 
     # Create mask model
@@ -270,7 +272,7 @@ def test_dqinit_refpix(instrument, exptype):
 
     # check if reference pixels are correct
     assert result.data.shape == (2, 20, 20)  # no pixels should be trimmed
-    assert result.amp33.shape == (2, 4096 ,128)
+    assert result.amp33.value.shape == (2, 4096 ,128)
     assert result.border_ref_pix_right.shape == (2, 20, 4)
     assert result.border_ref_pix_left.shape == (2, 20, 4)
     assert result.border_ref_pix_top.shape == (2, 4, 20)
diff --git a/romancal/flatfield/flat_field.py b/romancal/flatfield/flat_field.py
index 1fdcc1800..429632e0d 100644
--- a/romancal/flatfield/flat_field.py
+++ b/romancal/flatfield/flat_field.py
@@ -4,6 +4,8 @@
 
 import logging
 import numpy as np
+from astropy import units as u
+from roman_datamodels import units as ru
 
 from romancal.lib import dqflags
 
@@ -109,21 +111,21 @@ def apply_flat_field(science, flat):
     flat_data[np.where(flat_bad)] = 1.0
     # Now let's apply the correction to science data and error arrays.  Rely
     # on array broadcasting to handle the cubes
-    science.data /= flat_data
+    science.data = u.Quantity((science.data.value / flat_data), ru.electron / u.s, dtype=science.data.dtype)
 
     # Update the variances using BASELINE algorithm.  For guider data, it has
     # not gone through ramp fitting so there is no Poisson noise or readnoise
     flat_data_squared = flat_data**2
-    science.var_poisson /= flat_data_squared
-    science.var_rnoise /= flat_data_squared
-    try:
-        science.var_flat = science.data**2 / flat_data_squared * flat_err**2
-    except AttributeError:
-        science['var_flat'] = np.zeros(shape=science.data.shape,
-                                       dtype=np.float32)
-        science.var_flat = science.data**2 / flat_data_squared * flat_err**2
-    science.err = np.sqrt(science.var_poisson +
-                          science.var_rnoise + science.var_flat)
+    science.var_poisson = u.Quantity((science.var_poisson.value / flat_data_squared), ru.electron**2 / u.s**2,
+                                     dtype = science.var_poisson.dtype)
+    science.var_rnoise = u.Quantity((science.var_rnoise / flat_data_squared), ru.electron**2 / u.s**2,
+                                    dtype = science.var_rnoise.dtype)
+
+    science['var_flat'] = u.Quantity((science.data.value ** 2 / flat_data_squared * flat_err ** 2),
+                                  ru.electron ** 2 / u.s ** 2, dtype=science.data.value.dtype)
+
+    err_sqrt = np.sqrt(science.var_poisson.value + science.var_rnoise.value + science.var_flat)
+    science.err = u.Quantity(err_sqrt, ru.electron / u.s, dtype=err_sqrt.dtype)
 
     # Combine the science and flat DQ arrays
     science.dq = np.bitwise_or(science.dq, flat_dq)
diff --git a/romancal/flatfield/tests/test_flatfield.py b/romancal/flatfield/tests/test_flatfield.py
index 83774307c..ac715cfd1 100644
--- a/romancal/flatfield/tests/test_flatfield.py
+++ b/romancal/flatfield/tests/test_flatfield.py
@@ -6,9 +6,10 @@
 from roman_datamodels import stnode
 from roman_datamodels.datamodels import ImageModel, FlatRefModel
 from roman_datamodels.testing import utils as testutil
+from roman_datamodels import units as ru
 from romancal.flatfield import FlatFieldStep
 from astropy.time import Time
-
+from astropy import units as u
 
 @pytest.mark.parametrize(
     "instrument, exptype",
@@ -31,12 +32,17 @@ def test_flatfield_step_interface(instrument, exptype):
     wfi_image.meta.instrument.detector = 'WFI01'
     wfi_image.meta.instrument.optical_element = 'F158'
     wfi_image.meta.exposure.type = exptype
-    wfi_image.data = np.ones(shape, dtype=np.float32)
+    wfi_image.data = u.Quantity(np.ones(shape, dtype=np.float32),
+                                ru.electron / u.s, dtype=np.float32)
     wfi_image.dq = np.zeros(shape, dtype=np.uint32)
-    wfi_image.err = np.zeros(shape, dtype=np.float32)
-    wfi_image.var_poisson = np.zeros(shape, dtype=np.float32)
-    wfi_image.var_rnoise = np.zeros(shape, dtype=np.float32)
-    wfi_image.var_flat = np.zeros(shape, dtype=np.float32)
+    wfi_image.err = u.Quantity(np.zeros(shape, dtype=np.float32),
+                               ru.electron / u.s, dtype=np.float32)
+    wfi_image.var_poisson = u.Quantity(np.zeros(shape, dtype=np.float32),
+                                       ru.electron**2 / u.s**2, dtype=np.float32)
+    wfi_image.var_rnoise = u.Quantity(np.zeros(shape, dtype=np.float32),
+                                      ru.electron**2 / u.s**2, dtype=np.float32)
+    wfi_image.var_flat = u.Quantity(np.zeros(shape, dtype=np.float32),
+                                    ru.electron**2 / u.s**2, dtype=np.float32)
     wfi_image_model = ImageModel(wfi_image)
     flatref = stnode.FlatRef()
     meta = {}
diff --git a/romancal/jump/jump_step.py b/romancal/jump/jump_step.py
index 259b5329e..3ee1cdd81 100644
--- a/romancal/jump/jump_step.py
+++ b/romancal/jump/jump_step.py
@@ -58,10 +58,10 @@ def process(self, input):
             frames_per_group = meta.exposure.nframes
 
             # Modify the arrays for input into the 'common' jump (4D)
-            data = np.copy(r_data[np.newaxis, :])
+            data = np.copy(r_data[np.newaxis, :].value)
             gdq = r_gdq[np.newaxis, :]
             pdq = r_pdq[np.newaxis, :]
-            err = np.copy(r_err[np.newaxis, :])
+            err = np.copy(r_err[np.newaxis, :].value)
 
             tstart = time.time()
 
diff --git a/romancal/jump/tests/test_jump_step.py b/romancal/jump/tests/test_jump_step.py
index 71f003d3b..5da7195e9 100644
--- a/romancal/jump/tests/test_jump_step.py
+++ b/romancal/jump/tests/test_jump_step.py
@@ -7,12 +7,14 @@
 import pytest
 import numpy as np
 from astropy.time import Time
+from astropy import units as u
 
 import roman_datamodels.stnode as rds
 from roman_datamodels import datamodels as rdm
 from roman_datamodels.datamodels import GainRefModel
 from roman_datamodels.datamodels import ReadnoiseRefModel
 from roman_datamodels.testing import utils as testutil
+from roman_datamodels import units as ru
 
 from romancal.jump import JumpStep
 
@@ -49,7 +51,7 @@ def generate_wfi_reffiles(tmpdir_factory):
     gain_ref['meta'] = meta
     gain_ref['data'] = np.ones(shape, dtype=np.float32) * ingain
     gain_ref['dq'] = np.zeros(shape, dtype=np.uint16)
-    gain_ref['err'] = (np.random.random(shape) * 0.05).astype(np.float64)
+    gain_ref['err'] = (np.random.random(shape) * 0.05).astype(np.float32)
 
     gain_ref_model = GainRefModel(gain_ref)
     gain_ref_model.save(gainfile)
@@ -72,7 +74,7 @@ def generate_wfi_reffiles(tmpdir_factory):
     rn_ref['meta'] = meta
     rn_ref['data'] = np.ones(shape, dtype=np.float32)
     rn_ref['dq'] = np.zeros(shape, dtype=np.uint16)
-    rn_ref['err'] = (np.random.random(shape) * 0.05).astype(np.float64)
+    rn_ref['err'] = (np.random.random(shape) * 0.05).astype(np.float32)
 
     rn_ref_model = ReadnoiseRefModel(rn_ref)
     rn_ref_model.save(readnoisefile)
@@ -101,10 +103,10 @@ def _setup(ngroups=10, readnoise=10, nrows=20, ncols=20,
         dm_ramp.meta.instrument.name = 'WFI'
         dm_ramp.meta.instrument.optical_element = 'F158'
 
-        dm_ramp.data = data + 6.
+        dm_ramp.data = u.Quantity(data + 6., ru.DN, dtype=data.dtype)
         dm_ramp.pixeldq = pixdq
         dm_ramp.groupdq = gdq
-        dm_ramp.err = err
+        dm_ramp.err = u.Quantity(err, ru.DN, dtype=err.dtype)
 
         dm_ramp.meta.exposure.type = 'WFI_IMAGE'
         dm_ramp.meta.exposure.group_time = deltatime
@@ -142,7 +144,7 @@ def test_one_CR(generate_wfi_reffiles, max_cores, setup_inputs):
                           deltatime=grouptime)
 
     for i in range(ngroups):
-        model1.data[i, :, :] = deltaDN * i
+        model1.data.value[i, :, :] = deltaDN * i
     first_CR_group_locs = [x for x in range(1, 7) if x % 5 == 0]
 
     CR_locs = [x for x in range(xsize * ysize) if x % CR_fraction == 0]
@@ -153,8 +155,8 @@ def test_one_CR(generate_wfi_reffiles, max_cores, setup_inputs):
     # Add CRs to the SCI data
     for i in range(len(CR_x_locs)):
         CR_group = next(CR_pool)
-        model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]] = \
-            model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]] + 500.
+        model1.data.value[CR_group:, CR_y_locs[i], CR_x_locs[i]] = \
+            model1.data.value[CR_group:, CR_y_locs[i], CR_x_locs[i]] + 500.
 
     out_model = JumpStep.call(model1, override_gain=override_gain,
                               override_readnoise=override_readnoise,
@@ -190,7 +192,7 @@ def test_two_CRs(generate_wfi_reffiles, max_cores, setup_inputs):
                           deltatime=grouptime)
 
     for i in range(ngroups):
-        model1.data[i, :, :] = deltaDN * i
+        model1.data.value[i, :, :] = deltaDN * i
 
     first_CR_group_locs = [x for x in range(1, 7) if x % 5 == 0]
     CR_locs = [x for x in range(xsize * ysize) if x % CR_fraction == 0]
@@ -201,10 +203,10 @@ def test_two_CRs(generate_wfi_reffiles, max_cores, setup_inputs):
     for i in range(len(CR_x_locs)):
         CR_group = next(CR_pool)
 
-        model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]] = \
-            model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]] + 500
-        model1.data[CR_group + 8:, CR_y_locs[i], CR_x_locs[i]] = \
-            model1.data[CR_group + 8:, CR_y_locs[i], CR_x_locs[i]] + 700
+        model1.data.value[CR_group:, CR_y_locs[i], CR_x_locs[i]] = \
+            model1.data.value[CR_group:, CR_y_locs[i], CR_x_locs[i]] + 500
+        model1.data.value[CR_group + 8:, CR_y_locs[i], CR_x_locs[i]] = \
+            model1.data.value[CR_group + 8:, CR_y_locs[i], CR_x_locs[i]] + 700
 
     out_model = JumpStep.call(model1, override_gain=override_gain,
                               override_readnoise=override_readnoise,
diff --git a/romancal/linearity/linearity_step.py b/romancal/linearity/linearity_step.py
index 13e7d0fe7..332fab87a 100644
--- a/romancal/linearity/linearity_step.py
+++ b/romancal/linearity/linearity_step.py
@@ -9,6 +9,7 @@
 from romancal.stpipe import RomanStep
 from romancal.lib import dqflags
 from stcal.linearity.linearity import linearity_correction
+from astropy import units as u
 
 __all__ = ["LinearityStep"]
 
@@ -56,11 +57,11 @@ def process(self, input):
             # Call linearity correction function in stcal
             # The third return value is the procesed zero frame which
             # Roman does not use.
-            new_data, new_pdq, _ = linearity_correction(output_model.data,
+            new_data, new_pdq, _ = linearity_correction(output_model.data.value,
                                                      gdq, pdq, lin_coeffs,
                                                      lin_dq, dqflags.pixel)
 
-            output_model.data = new_data[0, :, :, :]
+            output_model.data = u.Quantity(new_data[0, :, :, :], u.DN, dtype=new_data.dtype)
             output_model.pixeldq = new_pdq
 
             # Close the reference file and update the step status
diff --git a/romancal/ramp_fitting/ramp_fit_step.py b/romancal/ramp_fitting/ramp_fit_step.py
index dc628cb8d..27436f1f5 100644
--- a/romancal/ramp_fitting/ramp_fit_step.py
+++ b/romancal/ramp_fitting/ramp_fit_step.py
@@ -8,6 +8,8 @@
 from roman_datamodels import datamodels as rdd
 from roman_datamodels import stnode as rds
 from roman_datamodels.testing import utils as testutil
+from roman_datamodels import units as ru
+from astropy import units as u
 
 from stcal.ramp_fitting import ramp_fit
 
@@ -42,15 +44,15 @@ def create_optional_results_model(input_model, opt_info):
     crmag.dtype = np.float32
 
     inst = {'meta': meta,
-            'slope': np.squeeze(slope),
-            'sigslope': np.squeeze(sigslope),
-            'var_poisson': np.squeeze(var_poisson),
-            'var_rnoise': np.squeeze(var_rnoise),
-            'yint': np.squeeze(yint),
-            'sigyint': np.squeeze(sigyint),
-            'pedestal': np.squeeze(pedestal),
+            'slope': u.Quantity(np.squeeze(slope), ru.electron / u.s, dtype=slope.dtype),
+            'sigslope': u.Quantity(np.squeeze(sigslope), ru.electron / u.s, dtype=sigslope.dtype),
+            'var_poisson': u.Quantity(np.squeeze(var_poisson), ru.electron**2 / u.s**2, dtype=var_poisson.dtype),
+            'var_rnoise': u.Quantity(np.squeeze(var_rnoise), ru.electron**2 / u.s**2, dtype=var_rnoise.dtype),
+            'yint': u.Quantity(np.squeeze(yint), ru.electron, dtype=yint.dtype),
+            'sigyint': u.Quantity(np.squeeze(sigyint), ru.electron, dtype=sigyint.dtype),
+            'pedestal': u.Quantity(np.squeeze(pedestal), ru.electron, dtype=pedestal.dtype),
             'weights': np.squeeze(weights),
-            'crmag': crmag
+            'crmag': u.Quantity(crmag, ru.electron, dtype=pedestal.dtype),
             }
 
     out_node = rds.RampFitOutput(inst)
@@ -80,6 +82,11 @@ def create_image_model(input_model, image_info):
     """
     data, dq, var_poisson, var_rnoise, err = image_info
 
+    data = u.Quantity(data, ru.electron / u.s, dtype=data.dtype)
+    var_poisson = u.Quantity(var_poisson, ru.electron**2 / u.s**2, dtype=var_poisson.dtype)
+    var_rnoise = u.Quantity(var_rnoise, ru.electron**2 / u.s**2, dtype=var_rnoise.dtype)
+    err = u.Quantity(err, ru.electron / u.s, dtype=err.dtype)
+
     # Create output datamodel
     # ... and add all keys from input
     meta = {}
@@ -142,7 +149,6 @@ def process(self, input):
             readnoise_filename = self.get_reference_file(input_model, 'readnoise')
             gain_filename = self.get_reference_file(input_model, 'gain')
             input_model.data = input_model.data[np.newaxis, :]
-            input_model.data.dtype=np.float32
             input_model.groupdq = input_model.groupdq[np.newaxis, :]
             input_model.err = input_model.err[np.newaxis, :]
 
diff --git a/romancal/ramp_fitting/tests/test_ramp_fit.py b/romancal/ramp_fitting/tests/test_ramp_fit.py
index 86d4549e9..f059f92e3 100644
--- a/romancal/ramp_fitting/tests/test_ramp_fit.py
+++ b/romancal/ramp_fitting/tests/test_ramp_fit.py
@@ -2,9 +2,11 @@
 import numpy as np
 import os
 from astropy.time import Time
+from astropy import units as u
 
 from roman_datamodels.datamodels import RampModel, GainRefModel, ReadnoiseRefModel, ImageModel
 from roman_datamodels.testing import utils as testutil
+from roman_datamodels import units as ru
 
 from romancal.ramp_fitting import RampFitStep
 from romancal.lib import dqflags
@@ -29,12 +31,10 @@ def generate_ramp_model(shape, deltatime=1):
     gdq = np.zeros(shape=shape, dtype=np.uint8)
 
     dm_ramp = testutil.mk_ramp(shape)
-    dm_ramp.data = data
+    dm_ramp.data = u.Quantity(data, ru.DN, dtype=np.float32)
     dm_ramp.pixeldq = pixdq
     dm_ramp.groupdq = gdq
-    dm_ramp.err = err
-
-    #dm_ramp.meta['photometry'] = testutil.mk_photometry()
+    dm_ramp.err = u.Quantity(err, ru.DN, dtype=np.float32)
 
     dm_ramp.meta.exposure.frame_time = deltatime
     dm_ramp.meta.exposure.ngroups = shape[0]
@@ -57,7 +57,7 @@ def generate_wfi_reffiles(shape, ingain = 6):
 
     gain_ref['data'] = (np.random.random(shape) * 0.5).astype(np.float32) * ingain
     gain_ref['dq'] = np.zeros(shape, dtype=np.uint16)
-    gain_ref['err'] = (np.random.random(shape) * 0.05).astype(np.float64)
+    gain_ref['err'] = (np.random.random(shape) * 0.05).astype(np.float32)
 
     gain_ref_model = GainRefModel(gain_ref)
 
@@ -95,7 +95,7 @@ def test_one_group_small_buffer_fit_ols(max_cores):
 
     model1 = generate_ramp_model(shape, deltatime)
 
-    model1.data[0, 15, 10] = 10.0  # add single CR
+    model1.data.value[0, 15, 10] = 10.0  # add single CR
 
     out_model = \
         RampFitStep.call(model1, override_gain=override_gain,
@@ -105,7 +105,7 @@ def test_one_group_small_buffer_fit_ols(max_cores):
     data = out_model.data
 
     # Index changes due to trimming of reference pixels
-    np.testing.assert_allclose(data[11, 6], 10.0, 1e-6)
+    np.testing.assert_allclose(data.value[11, 6], 10.0, 1e-6)
 
 @pytest.mark.skipif(
     os.environ.get("CI") == "true",
@@ -178,10 +178,10 @@ def test_saturated_ramp_fit(max_cores):
                          maximum_cores=max_cores)
 
     # Test data and error arrays are zeroed out
-    np.testing.assert_array_equal(out_model.data, 0)
-    np.testing.assert_array_equal(out_model.err, 0)
-    np.testing.assert_array_equal(out_model.var_poisson, 0)
-    np.testing.assert_array_equal(out_model.var_rnoise, 0)
+    np.testing.assert_array_equal(out_model.data.value, 0)
+    np.testing.assert_array_equal(out_model.err.value, 0)
+    np.testing.assert_array_equal(out_model.var_poisson.value, 0)
+    np.testing.assert_array_equal(out_model.var_rnoise.value, 0)
 
     # Test that all pixels are flagged saturated
     assert np.all(np.bitwise_and(out_model.dq, SATURATED) == SATURATED)
diff --git a/romancal/saturation/saturation.py b/romancal/saturation/saturation.py
index 4fb9b95b1..b7ae5daaf 100644
--- a/romancal/saturation/saturation.py
+++ b/romancal/saturation/saturation.py
@@ -34,7 +34,7 @@ def flag_saturation(input_model, ref_model):
         the GROUPDQ array
     """
 
-    data = input_model.data[np.newaxis, :]
+    data = input_model.data[np.newaxis, :].value
 
     # Create the output model as a copy of the input
     output_model = input_model
diff --git a/romancal/saturation/tests/test_saturation.py b/romancal/saturation/tests/test_saturation.py
index 16ed19106..5c62bd2c3 100644
--- a/romancal/saturation/tests/test_saturation.py
+++ b/romancal/saturation/tests/test_saturation.py
@@ -21,23 +21,23 @@ def test_basic_saturation_flagging(setup_wfi_datamodels):
     nrows = 20
     ncols = 20
     satvalue = 60000
-    data, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
+    ramp, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
 
     # Add ramp values up to the saturation limit
-    data.data[0, 5, 5] = 0
-    data.data[1, 5, 5] = 20000
-    data.data[2, 5, 5] = 40000
-    data.data[3, 5, 5] = 60000   # Signal reaches saturation limit
-    data.data[4, 5, 5] = 62000
+    ramp.data.value[0, 5, 5] = 0
+    ramp.data.value[1, 5, 5] = 20000
+    ramp.data.value[2, 5, 5] = 40000
+    ramp.data.value[3, 5, 5] = 60000   # Signal reaches saturation limit
+    ramp.data.value[4, 5, 5] = 62000
 
     # Set saturation value in the saturation model
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap)
+    output = flag_saturation(ramp, satmap)
 
     # Make sure that groups with signal > saturation limit get flagged
-    satindex = np.argmax(output.data[:, 5, 5] == satvalue)
+    satindex = np.argmax(output.data.value[:, 5, 5] == satvalue)
     assert np.all(output.groupdq[satindex:, 5, 5] == dqflags.group['SATURATED'])
 
 def test_ad_floor_flagging(setup_wfi_datamodels):
@@ -50,14 +50,14 @@ def test_ad_floor_flagging(setup_wfi_datamodels):
     ncols = 20
     satvalue = 60000
 
-    data, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
+    ramp, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
 
     # Add ramp values up to the saturation limit
-    data.data[0, 5, 5] = 0  # Signal at bottom rail - low saturation
-    data.data[1, 5, 5] = 0  # Signal at bottom rail - low saturation
-    data.data[2, 5, 5] = 20
-    data.data[3, 5, 5] = 40
-    data.data[4, 5, 5] = 60
+    ramp.data.value[0, 5, 5] = 0  # Signal at bottom rail - low saturation
+    ramp.data.value[1, 5, 5] = 0  # Signal at bottom rail - low saturation
+    ramp.data.value[2, 5, 5] = 20
+    ramp.data.value[3, 5, 5] = 40
+    ramp.data.value[4, 5, 5] = 60
 
     # frames that should be flagged as saturation (low)
     satindxs = [0, 1]
@@ -66,7 +66,7 @@ def test_ad_floor_flagging(setup_wfi_datamodels):
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap)
+    output = flag_saturation(ramp, satmap)
 
     # Check if the right frames are flagged as saturated
     assert np.all(output.groupdq[satindxs, 5, 5]
@@ -83,14 +83,14 @@ def test_ad_floor_and_saturation_flagging(setup_wfi_datamodels):
     ncols = 20
     satvalue = 60000
 
-    data, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
+    ramp, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
 
     # Add ramp values up to the saturation limit
-    data.data[0, 5, 5] = 0  # Signal at bottom rail - low saturation
-    data.data[1, 5, 5] = 0  # Signal at bottom rail - low saturation
-    data.data[2, 5, 5] = 20
-    data.data[3, 5, 5] = 40
-    data.data[4, 5, 5] = 61000  # Signal above the saturation threshold
+    ramp.data.value[0, 5, 5] = 0  # Signal at bottom rail - low saturation
+    ramp.data.value[1, 5, 5] = 0  # Signal at bottom rail - low saturation
+    ramp.data.value[2, 5, 5] = 20
+    ramp.data.value[3, 5, 5] = 40
+    ramp.data.value[4, 5, 5] = 61000  # Signal above the saturation threshold
 
     # frames that should be flagged as ad_floor
     floorindxs = [0, 1]
@@ -101,7 +101,7 @@ def test_ad_floor_and_saturation_flagging(setup_wfi_datamodels):
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap)
+    output = flag_saturation(ramp, satmap)
 
     # Check if the right frames are flagged as ad_floor
     assert np.all(output.groupdq[floorindxs, 5, 5]
@@ -121,23 +121,23 @@ def test_signal_fluctuation_flagging(setup_wfi_datamodels):
     ncols = 20
     satvalue = 60000
 
-    data, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
+    ramp, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
 
     # Add ramp values up to the saturation limit
-    data.data[0, 5, 5] = 10
-    data.data[1, 5, 5] = 20000
-    data.data[2, 5, 5] = 40000
-    data.data[3, 5, 5] = 60000   # Signal reaches saturation limit
-    data.data[4, 5, 5] = 40000   # Signal drops below saturation limit
+    ramp.data.value[0, 5, 5] = 10
+    ramp.data.value[1, 5, 5] = 20000
+    ramp.data.value[2, 5, 5] = 40000
+    ramp.data.value[3, 5, 5] = 60000   # Signal reaches saturation limit
+    ramp.data.value[4, 5, 5] = 40000   # Signal drops below saturation limit
 
     # Set saturation value in the saturation model
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap)
+    output = flag_saturation(ramp, satmap)
 
     # Make sure that all groups after first saturated group are flagged
-    satindex = np.argmax(output.data[:, 5, 5] == satvalue)
+    satindex = np.argmax(output.data.value[:, 5, 5] == satvalue)
     assert np.all(output.groupdq[satindex:, 5, 5] == dqflags.group['SATURATED'])
 
 
@@ -150,20 +150,20 @@ def test_all_groups_saturated(setup_wfi_datamodels):
     ncols = 20
     satvalue = 60000
 
-    data, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
+    ramp, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
 
     # Add ramp values at or above saturation limit
-    data.data[0, 5, 5] = 60000
-    data.data[1, 5, 5] = 62000
-    data.data[2, 5, 5] = 62000
-    data.data[3, 5, 5] = 60000
-    data.data[4, 5, 5] = 62000
+    ramp.data.value[0, 5, 5] = 60000
+    ramp.data.value[1, 5, 5] = 62000
+    ramp.data.value[2, 5, 5] = 62000
+    ramp.data.value[3, 5, 5] = 60000
+    ramp.data.value[4, 5, 5] = 62000
 
     # Set saturation value in the saturation model
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap)
+    output = flag_saturation(ramp, satmap)
 
     # Make sure all groups are flagged
     assert np.all(output.groupdq[:, 5, 5] == dqflags.group['SATURATED'])
@@ -178,14 +178,14 @@ def test_dq_propagation(setup_wfi_datamodels):
     dqval1 = 5
     dqval2 = 10
 
-    data, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
+    ramp, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
 
     # Add DQ values to the data and reference file
-    data.pixeldq[5, 5] = dqval1
+    ramp.pixeldq[5, 5] = dqval1
     satmap.dq[5, 5] = dqval2
 
     # Run the pipeline
-    output = flag_saturation(data, satmap)
+    output = flag_saturation(ramp, satmap)
 
     # Make sure DQ values from data and reference file are added in the output
     assert output.pixeldq[5, 5] == dqval1 + dqval2
@@ -201,24 +201,24 @@ def test_no_sat_check(setup_wfi_datamodels):
     ncols = 20
     satvalue = 60000
 
-    data, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
+    ramp, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
 
     # Add ramp values up to the saturation limit
-    data.data[0, 5, 5] = 10
-    data.data[1, 5, 5] = 20000
-    data.data[2, 5, 5] = 40000
-    data.data[3, 5, 5] = 60000
-    data.data[4, 5, 5] = 62000   # Signal reaches saturation limit
+    ramp.data.value[0, 5, 5] = 10
+    ramp.data.value[1, 5, 5] = 20000
+    ramp.data.value[2, 5, 5] = 40000
+    ramp.data.value[3, 5, 5] = 60000
+    ramp.data.value[4, 5, 5] = 62000   # Signal reaches saturation limit
 
     # Set saturation value in the saturation model & DQ value for NO_SAT_CHECK
     satmap.data[5, 5] = satvalue
     satmap.dq[5, 5] = dqflags.pixel['NO_SAT_CHECK']
 
     # Also set an existing DQ flag in input science data
-    data.pixeldq[5, 5] = dqflags.pixel['DO_NOT_USE']
+    ramp.pixeldq[5, 5] = dqflags.pixel['DO_NOT_USE']
 
     # Run the pipeline
-    output = flag_saturation(data, satmap)
+    output = flag_saturation(ramp, satmap)
 
     # Make sure output GROUPDQ does not get flagged as saturated
     # Make sure PIXELDQ is set to NO_SAT_CHECK and original flag
@@ -239,20 +239,20 @@ def test_nans_in_mask(setup_wfi_datamodels):
     nrows = 20
     ncols = 20
 
-    data, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
+    ramp, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
 
     # Add ramp values up to the saturation limit
-    data.data[0, 5, 5] = 10
-    data.data[1, 5, 5] = 20000
-    data.data[2, 5, 5] = 40000
-    data.data[3, 5, 5] = 60000
-    data.data[4, 5, 5] = 62000
+    ramp.data.value[0, 5, 5] = 10
+    ramp.data.value[1, 5, 5] = 20000
+    ramp.data.value[2, 5, 5] = 40000
+    ramp.data.value[3, 5, 5] = 60000
+    ramp.data.value[4, 5, 5] = 62000
 
     # Set saturation value for pixel to NaN
     satmap.data[5, 5] = np.nan
 
     # Run the pipeline
-    output = flag_saturation(data, satmap)
+    output = flag_saturation(ramp, satmap)
 
     # Check that output GROUPDQ is not flagged as saturated
     assert np.all(output.groupdq[:, 5, 5] != dqflags.group['SATURATED'])