diff --git a/CHANGES.rst b/CHANGES.rst
index 21de11519..2cfbbe992 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -12,6 +12,12 @@ Documentation
 - Add documentation for DNS build 0.5, e.g. reference array trimming [#457]
 
 
+photom
+------
+
+ - Added photom correction step and unit tests. [#469]
+
+
 0.6.0 (2022-03-02)
 ==================
 
diff --git a/romancal/lib/suffix.py b/romancal/lib/suffix.py
index 30b9722d2..4aa677f7a 100644
--- a/romancal/lib/suffix.py
+++ b/romancal/lib/suffix.py
@@ -67,6 +67,7 @@
     'assign_wcs',
     'linearity',
     'rampfitstep',
+    'photomstep',
     'pipeline',
     'dq_init',
     'linearitystep',
diff --git a/romancal/photom/__init__.py b/romancal/photom/__init__.py
new file mode 100644
index 000000000..59fa0f31b
--- /dev/null
+++ b/romancal/photom/__init__.py
@@ -0,0 +1,3 @@
+from .photom_step import PhotomStep
+
+__all__ = ['PhotomStep']
diff --git a/romancal/photom/photom.py b/romancal/photom/photom.py
new file mode 100644
index 000000000..06385fc3b
--- /dev/null
+++ b/romancal/photom/photom.py
@@ -0,0 +1,110 @@
+import logging
+import warnings
+from astropy import units as u
+
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
+
+def photom_io(input_model, photom_metadata):
+    """
+    Combine photometric scalar conversion factors and add to the science metadata.
+
+    Parameters
+    ----------
+    input_model : Roman level 2 image datamodel
+        Input Roman datamodel
+
+    photom_metadata : dict
+        Set of matched photom meta data keys
+
+    Returns
+    -------
+
+    """
+    # Get the scalar conversion factor.
+    conversion = photom_metadata['photmjsr']  # unit is MJy / sr
+
+    # Store the conversion factor in the meta data
+    log.info(f'photmjsr value: {conversion:.6g}')
+    input_model.meta.photometry.conversion_megajanskys = conversion
+    input_model.meta.photometry.conversion_microjanskys = conversion.to(
+        u.microjansky / u.arcsecond ** 2)
+
+    # Get the scalar conversion uncertainty factor
+    uncertainty_conv = photom_metadata['uncertainty']
+
+    # Store the uncertainty conversion factor in the meta data
+    log.info(f'uncertainty value: {conversion:.6g}')
+    input_model.meta.photometry.conversion_megajanskys_uncertainty = uncertainty_conv
+    input_model.meta.photometry.conversion_microjanskys_uncertainty = uncertainty_conv.to(
+        u.microjansky / u.arcsecond ** 2)
+
+    # Return updated input model
+    return input_model
+
+
+def save_area_info(input_model, photom_parameters):
+    """
+    Read the pixel area value in the photom parameters, then convert and
+    copy them to the metadata of the input datamodel.
+
+    Parameters
+    ----------
+    input_model : Roman level 2 image datamodel
+        Input Roman datamodel
+
+    photom_parameters : dict
+        Photom parameter object
+    """
+
+    # Load the average pixel area values from the photom reference file header
+    area_ster = photom_parameters["pixelareasr"]
+    area_a2 = photom_parameters["pixelareasr"].to(u.arcsecond**2)
+
+    # Copy the pixel area values to the input model
+    log.debug(f'pixelarea_steradians = {area_ster}, pixelarea_arcsecsq = {area_a2}')
+    input_model.meta.photometry.pixelarea_arcsecsq = area_a2
+    input_model.meta.photometry.pixelarea_steradians = area_ster
+
+    # Return updated input model
+    return input_model
+
+
+def apply_photom(input_model, photom):
+    """
+    Retrieve the conversion factor from the photom reference datamodel
+    that is appropriate to the instrument mode. The scalar factor is written
+    to the photmjsr and uncertainty keywords in the model.
+
+    For WFI, matching is based on optical_element.
+
+    Parameters
+    ----------
+    input_model : Roman level 2 image datamodel
+        Input Roman datamodel
+
+    photom : Roman Photom reference datamodel
+        Photom Roman datamodel
+
+    Returns
+    -------
+    output_model : Roman level 2 image datamodel
+        Output Roman datamodel with the flux calibration keywords updated
+    """
+    # Obtain photom parameters for the selected optical element
+    try:
+        photom_parameters = photom.phot_table[input_model.meta.instrument.optical_element.upper()]
+    except KeyError:
+        warnings.warn(f'No matching photom parameters for '
+                      f'{input_model.meta.instrument.optical_element}')
+        return input_model
+
+    # Copy pixel area information to output datamodel
+    output_model = save_area_info(input_model, photom_parameters)
+
+    # Copy conversions to output model
+    output_model = photom_io(output_model, photom_parameters)
+
+    # Return updated output model
+    return output_model
diff --git a/romancal/photom/photom_step.py b/romancal/photom/photom_step.py
new file mode 100644
index 000000000..f89bf3252
--- /dev/null
+++ b/romancal/photom/photom_step.py
@@ -0,0 +1,62 @@
+#! /usr/bin/env python
+
+from romancal.stpipe import RomanStep
+from romancal.photom import photom
+import roman_datamodels as rdm
+
+__all__ = ["PhotomStep"]
+
+
+class PhotomStep(RomanStep):
+    """
+    PhotomStep: Module for loading photometric conversion information from
+        reference files and attaching to the input science data model
+    """
+
+    reference_file_types = ['photom']
+
+    def process(self, input):
+        """Perform the photometric calibration step
+
+        Parameters
+        ----------
+        input : Roman level 2 image datamodel (wfi_image-1.x.x)
+            input roman datamodel
+
+        Returns
+        -------
+        output_model : Roman level 2 image datamodel (wfi_image-1.x.x)
+            output roman datamodel
+        """
+
+        # Open the input data model
+        with rdm.open(input) as input_model:
+
+            # Get reference file
+            reffile = self.get_reference_file(input_model, "photom")
+
+            # Open the relevant photom reference file as a datamodel
+            if reffile is not None:
+                # If there is a reference file, perform photom application
+                photom_model = rdm.open(reffile)
+                self.log.debug(f'Using PHOTOM ref file: {reffile}')
+
+                # Do the correction
+                output_model = photom.apply_photom(input_model, photom_model)
+                output_model.meta.cal_step.photom = 'COMPLETE'
+
+            else:
+                # Skip Photom step if no photom file
+                self.log.warning('No PHOTOM reference file found')
+                self.log.warning('Photom step will be skipped')
+                input_model.meta.cal_step.photom = 'SKIPPED'
+                return input_model
+
+        # Close the input and reference files
+        input_model.close()
+        try:
+            photom_model.close()
+        except AttributeError:
+            pass
+
+        return output_model
diff --git a/romancal/photom/tests/__init__.py b/romancal/photom/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/romancal/photom/tests/test_photom.py b/romancal/photom/tests/test_photom.py
new file mode 100644
index 000000000..394488b0c
--- /dev/null
+++ b/romancal/photom/tests/test_photom.py
@@ -0,0 +1,211 @@
+import os
+import pytest
+import numpy as np
+import warnings
+from astropy import units as u
+
+from romancal.photom import photom, PhotomStep
+from roman_datamodels.datamodels import ImageModel, WfiImgPhotomRefModel
+from roman_datamodels.testing import utils as testutil
+
+
+def create_photom_wfi_image(min_r=3.1, delta=0.1):
+    """Create a photom table for WFI.
+
+    Parameters
+    ----------
+    min_r : float
+        Minimum value to assign when populating the photometry array.
+
+    max_r : float
+        Maximum value to assign when populating the photometry array.
+
+    Returns
+    -------
+    photom_model : Roman Photom reference datamodel
+        photom roman datamodel
+
+    """
+
+    # Keyword and eleents list for phot_table construction
+    optical_element = ["F062", "F087", "F106", "F129", "W146", "F158", "F184", "F213",
+                       "GRISM", "PRISM", "DARK"]
+    none_type_elements = ["GRISM", "PRISM", "DARK"]
+    keyword = ["photmjsr", "uncertainty", "pixelareasr"]
+    nrows = len(optical_element)
+
+    # Create sample photometry keyword values
+    photmjsr = np.linspace(min_r, min_r + (nrows - 1.) * delta, nrows) * u.megajansky / u.steradian
+    uncertainty = np.linspace(min_r/20.0, min_r/20.0 + (nrows - 1.) * delta/20.0, nrows) * \
+                  u.megajansky / u.steradian
+
+    # Create sample area keyword values
+    area_ster = 2.31307642258977E-14 * u.steradian
+    pixelareasr = np.ones(nrows, dtype=np.float32) * area_ster
+
+    # Bundle values into a list
+    values = list(zip(photmjsr, uncertainty, pixelareasr))
+
+    #Create dictionary containing all values
+    reftab = {}
+    for element_idx, element in enumerate(optical_element):
+        key_dict = {}
+        for key_idx, key in enumerate(keyword):
+            key_dict[key] = values[element_idx][key_idx]
+        # GRISM, PRISM, and DARK optical elements shuld have their
+        # photomeetry keywords set to None
+        if element in none_type_elements:
+            key_dict["photmjsr"] = None
+            key_dict["uncertainty"] = None
+        reftab[element] = key_dict
+
+    # Create default datamodel
+    photom_model = testutil.mk_wfi_img_photom()
+
+    # Copy values above into defautl datamodel
+    photom_model.phot_table=reftab
+
+    return photom_model
+
+
+def test_no_photom_match():
+    """Test apply_photom warning for no match"""
+
+    # Create sample WFI Level 2 science datamodel
+    input_model = testutil.mk_level2_image()
+
+    # Create photom reference datamodel
+    photom_model = create_photom_wfi_image(min_r=3.1, delta=0.1)
+
+    # Remove key for failed test (that won't fail validation)
+    photom_model.phot_table.pop("W146")
+
+    # Select optical element
+    input_model.meta.instrument.optical_element = "W146"
+
+    # Set bad values which would be overwritten by apply_photom
+    input_model.meta.photometry.pixelarea_steradians = -1.0 * u.sr
+    input_model.meta.photometry.conversion_megajanskys = -1.0 * u.megajansky / u.steradian
+    input_model.meta.photometry.conversion_microjanskys_uncertainty = \
+        -1.0 * u.microjansky / u.arcsecond ** 2
+
+    with warnings.catch_warnings(record=True) as caught:
+        # Look for now non existent W146 optical element
+        output_model = photom.apply_photom(input_model, photom_model)
+
+        # Assert warning key matches that of the input file
+        assert str(caught[0].message).split()[-1] == input_model.meta.instrument.optical_element
+
+        # Assert that photom elements are not updated
+        assert output_model.meta.photometry.pixelarea_steradians == -1.0 * u.sr
+        assert output_model.meta.photometry.conversion_megajanskys == \
+               -1.0 * u.megajansky / u.steradian
+        assert output_model.meta.photometry.conversion_microjanskys_uncertainty == \
+               -1.0 * u.microjansky / u.arcsecond ** 2
+
+
+def test_apply_photom1():
+    """Test apply_photom applies correct metadata"""
+
+    # Create sample WFI Level 2 science datamodel
+    input_model = testutil.mk_level2_image()
+
+    # Create photom reference datamodel
+    photom_model = create_photom_wfi_image(min_r=3.1, delta=0.1)
+
+    # Select optical element
+    input_model.meta.instrument.optical_element = "W146"
+
+    # Apply photom correction for optical element W146
+    output_model = photom.apply_photom(input_model, photom_model)
+
+    # Set reference photometry
+    area_ster = 2.31307642258977E-14 * u.steradian
+    area_a2 = 0.000984102303070964 * u.arcsecond * u.arcsecond
+
+    # Tests for pixel areas
+    assert(np.isclose(output_model.meta.photometry.pixelarea_steradians.value,
+                        area_ster.value, atol=1.e-7))
+    assert output_model.meta.photometry.pixelarea_steradians.unit == area_ster.unit
+    assert(np.isclose(output_model.meta.photometry.pixelarea_arcsecsq.value,
+                        area_a2.value, atol=1.e-7))
+    assert output_model.meta.photometry.pixelarea_arcsecsq.unit == area_a2.unit
+
+    # Set reference photometry
+    phot_ster = 3.5 * u.megajansky / u.steradian
+    phot_a2 = phot_ster.to(u.microjansky / u.arcsecond ** 2)
+
+    # Tests for photometry
+    assert (np.isclose(output_model.meta.photometry.conversion_megajanskys.value,
+                       phot_ster.value, atol=1.e-7))
+    assert output_model.meta.photometry.conversion_megajanskys.unit == phot_ster.unit
+    assert (np.isclose(output_model.meta.photometry.conversion_microjanskys.value,
+                       phot_a2.value, atol=1.e-7))
+    assert output_model.meta.photometry.conversion_microjanskys.unit == phot_a2.unit
+
+    # Set reference photometric uncertainty
+    muphot_ster = 0.175 * u.megajansky / u.steradian
+    muphot_a2 = muphot_ster.to(u.microjansky / u.arcsecond ** 2)
+
+    # Tests for photometric uncertainty
+    assert (np.isclose(output_model.meta.photometry.conversion_megajanskys_uncertainty.value,
+                       muphot_ster.value, atol=1.e-7))
+    assert output_model.meta.photometry.conversion_megajanskys_uncertainty.unit == muphot_ster.unit
+    assert (np.isclose(output_model.meta.photometry.conversion_microjanskys_uncertainty.value,
+                       muphot_a2.value, atol=1.e-7))
+    assert output_model.meta.photometry.conversion_microjanskys_uncertainty.unit == muphot_a2.unit
+
+
+def test_apply_photom2():
+    """Test apply_photom does not change data values"""
+
+    # Create sample WFI Level 2 science datamodel
+    input_model = testutil.mk_level2_image()
+
+    # Create photom reference datamodel
+    photom_model = create_photom_wfi_image(min_r=3.1, delta=0.1)
+
+    # Apply photom correction
+    output_model = photom.apply_photom(input_model, photom_model)
+
+    # Select pixel for comparison
+    shape = input_model.data.shape
+    ix = shape[1] // 2
+    iy = shape[0] // 2
+
+    # Test that the data has not changed
+    assert (np.allclose(output_model.data[iy, ix], input_model.data[iy, ix], rtol=1.e-7))
+
+
+@pytest.mark.parametrize(
+    "instrument, exptype",
+    [
+        ("WFI", "WFI_IMAGE"),
+    ]
+)
+@pytest.mark.skipif(
+    os.environ.get("CI") == "true",
+    reason="Roman CRDS servers are not currently available outside the internal network"
+)
+def test_photom_step_interface(instrument, exptype):
+    """Test that the basic inferface works for data requiring a photom reffile"""
+
+    # Create a small area for the file
+    shape = (20, 20)
+
+    # Create input model
+    wfi_image = testutil.mk_level2_image(shape=shape)
+    wfi_image_model = ImageModel(wfi_image)
+
+    # Create photom model
+    photom = testutil.mk_wfi_img_photom()
+    photom_model = WfiImgPhotomRefModel(photom)
+
+    photom_model
+
+    # Run photom correction step
+    result = PhotomStep.call(wfi_image_model, override_photom=photom_model)
+
+    assert (result.data == wfi_image.data).all()
+    assert result.data.shape == shape
+    assert result.meta.cal_step.photom == 'COMPLETE'
diff --git a/romancal/pipeline/exposure_pipeline.py b/romancal/pipeline/exposure_pipeline.py
index 85d8d1822..bfe64a51e 100644
--- a/romancal/pipeline/exposure_pipeline.py
+++ b/romancal/pipeline/exposure_pipeline.py
@@ -11,6 +11,7 @@
 from romancal.jump import jump_step
 from romancal.dark_current import DarkCurrentStep
 from romancal.linearity import LinearityStep
+from romancal.photom import PhotomStep
 from romancal.ramp_fitting import ramp_fit_step
 from romancal.saturation import SaturationStep
 
@@ -44,7 +45,8 @@ class ExposurePipeline(RomanPipeline):
                  'jump': jump_step.JumpStep,
                  'rampfit': ramp_fit_step.RampFitStep,
                  'assign_wcs': AssignWcsStep,
-                 'flatfield': FlatFieldStep
+                 'flatfield': FlatFieldStep,
+                 'photom': PhotomStep,
                  }
 
     # start the actual processing
@@ -80,6 +82,7 @@ def process(self, input):
         else:
             log.info('Flat Field step is being SKIPPED')
             result.meta.cal_step.flat_field = 'SKIPPED'
+        result = self.photom(result)
 
         # setup output_file for saving
         self.setup_output(result)
diff --git a/romancal/step.py b/romancal/step.py
index 150961bd4..192a9bd47 100644
--- a/romancal/step.py
+++ b/romancal/step.py
@@ -7,6 +7,7 @@
 from .flatfield.flat_field_step import FlatFieldStep
 from .jump.jump_step import JumpStep
 from .linearity.linearity_step import LinearityStep
+from .photom.photom_step import PhotomStep
 from .ramp_fitting.ramp_fit_step import RampFitStep
 from .saturation.saturation_step import SaturationStep
 from .assign_wcs.assign_wcs_step import AssignWcsStep
@@ -18,6 +19,7 @@
     "FlatFieldStep",
     "JumpStep",
     "LinearityStep",
+    "PhotomStep",
     "RampFitStep",
     "SaturationStep",
     "AssignWcsStep",
diff --git a/romancal/stpipe/integration.py b/romancal/stpipe/integration.py
index 054593e53..70dff39ae 100644
--- a/romancal/stpipe/integration.py
+++ b/romancal/stpipe/integration.py
@@ -27,6 +27,7 @@ def get_steps():
         ("romancal.step.FlatFieldStep", None, False),
         ("romancal.step.JumpStep", None, False),
         ("romancal.step.LinearityStep", None, False),
+        ("romancal.step.PhotomStep", None, False),
         ("romancal.step.RampFitStep", None, False),
         ("romancal.step.SaturationStep", None, False),
         ("romancal.step.AssignWcsStep", None, False),