diff --git a/lasy/laser.py b/lasy/laser.py
index 1c240c8f..36273daa 100644
--- a/lasy/laser.py
+++ b/lasy/laser.py
@@ -4,6 +4,7 @@
 
 from lasy.utils.grid import Grid, time_axis_indx
 from lasy.utils.laser_utils import (
+    normalize_average_intensity,
     normalize_energy,
     normalize_peak_field_amplitude,
     normalize_peak_intensity,
@@ -147,7 +148,7 @@ def __init__(
 
     def normalize(self, value, kind="energy"):
         """
-        Normalize the pulse either to the energy, peak field amplitude or peak intensity.
+        Normalize the pulse either to the energy, peak field amplitude, peak intensity, or average intensity. The average intensity option operates on the envelope.
 
         Parameters
         ----------
@@ -155,7 +156,7 @@ def normalize(self, value, kind="energy"):
             Value to which to normalize the field property that is defined in ``kind``
         kind : string (optional)
             Distance by which the laser pulse should be propagated
-            Options: ``'energy``', ``'field'``, ``'intensity'`` (default is ``'energy'``)
+            Options: ``'energy``', ``'field'``, ``'intensity'``, ``'average_intensity'`` (default is ``'energy'``)
         """
         if kind == "energy":
             normalize_energy(self.dim, value, self.grid)
@@ -163,6 +164,8 @@ def normalize(self, value, kind="energy"):
             normalize_peak_field_amplitude(value, self.grid)
         elif kind == "intensity":
             normalize_peak_intensity(value, self.grid)
+        elif kind == "average_intensity":
+            normalize_average_intensity(value, self.grid)
         else:
             raise ValueError(f'kind "{kind}" not recognized')
 
diff --git a/lasy/utils/laser_utils.py b/lasy/utils/laser_utils.py
index 55be1d3c..f9de3aa7 100644
--- a/lasy/utils/laser_utils.py
+++ b/lasy/utils/laser_utils.py
@@ -118,7 +118,7 @@ def normalize_peak_intensity(peak_intensity, grid):
     Parameters
     ----------
     peak_intensity : scalar (W/m^2)
-        Peak field amplitude of the laser pulse after normalization.
+        Peak intensity of the laser pulse after normalization.
 
     grid : a Grid object
         Contains value of the laser envelope and metadata.
@@ -130,7 +130,31 @@ def normalize_peak_intensity(peak_intensity, grid):
         if input_peak_intensity == 0.0:
             print("Field is zero everywhere, normalization will be skipped")
         else:
-            grid.set_temporal_field(np.sqrt(peak_intensity / input_peak_intensity))
+            field *= np.sqrt(peak_intensity / input_peak_intensity)
+            grid.set_temporal_field(field)
+
+
+def normalize_average_intensity(average_intensity, grid):
+    """
+    Normalize energy of the laser pulse contained in grid.
+
+    Parameters
+    ----------
+    average_intensity : scalar (W/m^2)
+        Average intensity of the laser pulse envelope after normalization.
+
+    grid : a Grid object
+        Contains value of the laser envelope and metadata.
+    """
+    if average_intensity is not None:
+        field = grid.get_temporal_field()
+        intensity = np.abs(epsilon_0 * field**2 / 2 * c)
+        input_average_intensity = intensity.mean()
+        if input_average_intensity == 0.0:
+            print("Field is zero everywhere, normalization will be skipped")
+        else:
+            field *= np.sqrt(average_intensity / input_average_intensity)
+            grid.set_temporal_field(field)
 
 
 def get_full_field(laser, theta=0, slice=0, slice_axis="x", Nt=None):
diff --git a/tests/test_laser_utils.py b/tests/test_laser_utils.py
index 0fbe926a..5446622f 100644
--- a/tests/test_laser_utils.py
+++ b/tests/test_laser_utils.py
@@ -1,4 +1,5 @@
 import numpy as np
+from scipy.constants import c, epsilon_0
 
 from lasy.laser import Laser
 from lasy.profiles.gaussian_profile import GaussianProfile
@@ -48,5 +49,34 @@ def test_laser_analysis_utils():
         np.testing.assert_approx_equal(2 * tau_rms, laser.profile.tau, significant=3)
 
 
+def test_laser_normalization_utils():
+    """Test the different laser normalization utilities in both geometries."""
+    for dim in ["xyt", "rt"]:
+        laser = get_gaussian_laser(dim)
+
+        # Check energy normalization
+        laser.normalize(1, kind="energy")
+        energy = compute_laser_energy(dim, laser.grid)
+        np.testing.assert_approx_equal(1, energy, significant=10)
+
+        # Check peak field normalization
+        laser.normalize(1, kind="field")
+        field = laser.grid.get_temporal_field()
+        np.testing.assert_approx_equal(1, np.abs(field.max()), significant=10)
+
+        # Check peak intensity normalization
+        laser.normalize(1, kind="intensity")
+        field = laser.grid.get_temporal_field()
+        intensity = np.abs(epsilon_0 * field**2 / 2 * c)
+        np.testing.assert_approx_equal(1, intensity.max(), significant=10)
+
+        # Check average intensity normalization
+        laser.normalize(1, kind="average_intensity")
+        field = laser.grid.get_temporal_field()
+        intensity = np.abs(epsilon_0 * field**2 / 2 * c)
+        np.testing.assert_approx_equal(1, intensity.mean(), significant=10)
+
+
 if __name__ == "__main__":
     test_laser_analysis_utils()
+    test_laser_normalization_utils()