diff --git a/Exec/Make.Castro b/Exec/Make.Castro
index 2db7d6aa36..e31c452842 100644
--- a/Exec/Make.Castro
+++ b/Exec/Make.Castro
@@ -166,6 +166,10 @@ ifeq ($(USE_CXX_MODEL_PARSER), TRUE)
   DEFINES += -DCXX_MODEL_PARSER -DNPTS_MODEL=$(MAX_NPTS_MODEL) -DNUM_MODELS=$(NUM_MODELS)
 endif
 
+ifeq ($(USE_RNG_STATE_INIT), TRUE)
+  DEFINES += -DRNG_STATE_INIT
+endif
+
 # add / define any special physics we need
 
 ifeq ($(USE_MHD), TRUE)
diff --git a/Exec/science/xrb_layered/GNUmakefile b/Exec/science/xrb_layered/GNUmakefile
index 516b853f83..f029cb400a 100644
--- a/Exec/science/xrb_layered/GNUmakefile
+++ b/Exec/science/xrb_layered/GNUmakefile
@@ -12,6 +12,7 @@ USE_GRAV         = TRUE
 USE_REACT        = TRUE
 
 USE_CXX_MODEL_PARSER = TRUE
+USE_RNG_STATE_INIT = TRUE
 
 CASTRO_HOME = ../../..
 
diff --git a/Exec/science/xrb_layered/problem_initialize_state_data.H b/Exec/science/xrb_layered/problem_initialize_state_data.H
index 8cc253f6b7..3eb8db2ecd 100644
--- a/Exec/science/xrb_layered/problem_initialize_state_data.H
+++ b/Exec/science/xrb_layered/problem_initialize_state_data.H
@@ -8,7 +8,8 @@
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void problem_initialize_state_data (int i, int j, int k,
                                     Array4<Real> const& state,
-                                    const GeometryData& geomdata)
+                                    const GeometryData& geomdata,
+                                    amrex::RandomEngine const& engine)
 {
 
     const Real* dx = geomdata.CellSize();
@@ -78,7 +79,7 @@ void problem_initialize_state_data (int i, int j, int k,
 
             // mean of 1, std dev of 1.e-3
 
-            state(i,j,k,UTEMP) *= amrex::RandomNormal(1.0, 0.001);
+            state(i,j,k,UTEMP) *= amrex::RandomNormal(1.0, 0.001, engine);
             eos_state.T = state(i,j,k,UTEMP);
 
             // recompute the density with this new temperature
diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp
index 7ef73d48e4..3e0ad948ea 100644
--- a/Source/driver/Castro.cpp
+++ b/Source/driver/Castro.cpp
@@ -1105,12 +1105,21 @@ Castro::initData ()
           auto s = S_new[mfi].array();
           auto geomdata = geom.data();
 
+#ifdef RNG_STATE_INIT
+          amrex::ParallelForRNG(box,
+          [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept
+          {
+              // problem initialization
+              problem_initialize_state_data(i, j, k, s, geomdata, engine);
+          });
+#else
           amrex::ParallelFor(box,
           [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
           {
               // problem initialization
               problem_initialize_state_data(i, j, k, s, geomdata);
           });
+#endif
        }
 
 
diff --git a/Source/driver/Castro_advance.cpp b/Source/driver/Castro_advance.cpp
index aa3e068748..a7190c2fd7 100644
--- a/Source/driver/Castro_advance.cpp
+++ b/Source/driver/Castro_advance.cpp
@@ -428,6 +428,9 @@ Castro::initialize_advance(Real time, Real dt, int amr_iteration)
                 auto s = S_old[mfi].array();
                 auto geomdata = geom.data();
 
+#ifdef RNG_STATE_INIT
+                amrex::Error("drive initial convection not yet supported for random initialization");
+#else
                 amrex::ParallelFor(box,
                 [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
                 {
@@ -449,6 +452,8 @@ Castro::initialize_advance(Real time, Real dt, int amr_iteration)
                         (vx_orig * vx_orig + vy_orig * vy_orig + vz_orig * vz_orig);
 
                 });
+#endif
+
             }
 
         }
diff --git a/Source/driver/Castro_io.cpp b/Source/driver/Castro_io.cpp
index e3bb00a5d2..33fde99403 100644
--- a/Source/driver/Castro_io.cpp
+++ b/Source/driver/Castro_io.cpp
@@ -291,6 +291,9 @@ Castro::restart (Amr&     papa,
                auto s = S_new[mfi].array();
                auto geomdata = geom.data();
 
+#ifdef RNG_STATE_INIT
+               amrex::Error("Error: random initialization not yet supported with grown factor");
+#else
                amrex::ParallelFor(bx,
                [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
                {
@@ -298,6 +301,7 @@ Castro::restart (Amr&     papa,
                    // by a problem setup (defaults to an empty routine).
                    problem_initialize_state_data(i, j, k, s, geomdata);
                });
+#endif
            }
        }
     }