diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90
index a1b02c1f86..a82978f2cf 100644
--- a/src/dynamics/mpas/dyn_comp.F90
+++ b/src/dynamics/mpas/dyn_comp.F90
@@ -11,7 +11,7 @@ module dyn_comp
 use const_init,         only: cnst_init_default
 
 use cam_control_mod,    only: initial_run
-use cam_initfiles,      only: initial_file_get_id, topo_file_get_id
+use cam_initfiles,      only: initial_file_get_id, topo_file_get_id, pertlim
 
 use cam_grid_support,   only: cam_grid_id, &
                               cam_grid_get_latvals, cam_grid_get_lonvals
@@ -813,6 +813,11 @@ subroutine read_inidat(dyn_in)
 
    logical  :: readvar
 
+   integer                          :: rndm_seed_sz
+   integer, allocatable             :: rndm_seed(:)
+   real(r8)                         :: pertval
+   integer                          :: nc
+
    character(len=shr_kind_cx) :: str
 
    type(mpas_pool_type), pointer :: mesh_pool
@@ -1083,6 +1088,29 @@ subroutine read_inidat(dyn_in)
          call endrun(subname//': failed to read theta from initial file')
       end if
 
+      ! optionally introduce random perturbations to theta values
+      if (pertlim.ne.0.0_r8) then
+         if (masterproc) then
+            write(iulog,*) trim(subname), ': Adding random perturbation bounded', &
+               'by +/- ', pertlim, ' to initial theta field'
+         end if
+
+         call random_seed(size=rndm_seed_sz)
+         allocate(rndm_seed(rndm_seed_sz))
+
+         do nc = 1,nCellsSolve
+            rndm_seed = glob_ind(nc)
+            call random_seed(put=rndm_seed)
+            do kk = 1,plev
+               call random_number(pertval)
+               pertval = 2.0_r8*pertlim*(0.5_r8 - pertval)
+               theta(kk,nc) = theta(kk,nc)*(1.0_r8 + pertval)
+            end do
+         end do
+
+         deallocate(rndm_seed)
+      end if
+
       ! read rho
       call infld('rho', fh_ini, 'lev', 'nCells', 1, plev, 1, nCellsSolve, 1, 1, &
                  mpas3d, readvar, gridname='mpas_cell')