diff --git a/Source/Initialization/ERF_Metgrid_utils.H b/Source/Initialization/ERF_Metgrid_utils.H
index 452e4f567..00a67cac4 100644
--- a/Source/Initialization/ERF_Metgrid_utils.H
+++ b/Source/Initialization/ERF_Metgrid_utils.H
@@ -3,7 +3,6 @@
 
 #include <ERF.H>
 #include <ERF_EOS.H>
-#include <ERF_Constants.H>
 #include <ERF_Utils.H>
 #include <ERF_prob_common.H>
 #include <ERF_HSE_utils.H>
@@ -107,6 +106,7 @@ void
 calc_rho_p (const int& kmax,
             const int& flag_psfc,
             const amrex::Real& psfc,
+            const amrex::Real& grav,
             amrex::Real* Thetad_vec,
             amrex::Real* Thetam_vec,
             amrex::Real* Q_vec,
@@ -126,7 +126,7 @@ calc_rho_p (const int& kmax,
     } else {
         amrex::Real t_0 = 290.0; // WRF's model_config_rec%base_temp
         amrex::Real a   = 50.0; // WRF's model_config_rec%base_lapse
-        Psurf = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*CONST_GRAV*z_vec[0]/(a*R_d)), 0.5));
+        Psurf = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*grav*z_vec[0]/(a*R_d)), 0.5));
     }
 
     // Iterations for the first CC point that is 1/2 dz off the surface
@@ -135,7 +135,7 @@ calc_rho_p (const int& kmax,
         amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[0];
         Thetam_vec[0] = Thetad_vec[0]*qvf;
         for (int it=0; it<maxiter; it++) {
-            Pm_vec[0] = Psurf-half_dz*(Rhom_vec[0])*(1.0+Q_vec[0])*CONST_GRAV;
+            Pm_vec[0] = Psurf-half_dz*(Rhom_vec[0])*(1.0+Q_vec[0])*grav;
             if (Pm_vec[0] < 0.0) Pm_vec[0] = 0.0;
             Rhom_vec[0] = (p_0/(R_d*Thetam_vec[0]))*std::pow(Pm_vec[0]/p_0, iGamma);
         } // it
@@ -153,7 +153,7 @@ calc_rho_p (const int& kmax,
 
         // Establish known constant
         amrex::Real rho_tot_lo = Rhom_vec[k-1] * (1. + Q_vec[k-1]);
-        amrex::Real C = -Pm_vec[k-1] + 0.5*rho_tot_lo*CONST_GRAV*dz;
+        amrex::Real C = -Pm_vec[k-1] + 0.5*rho_tot_lo*grav*dz;
 
         // Initial guess and residual
         Pm_vec[k]   = Pm_vec[k-1];
@@ -162,11 +162,11 @@ calc_rho_p (const int& kmax,
                                             R_d/Cp_d,
                                             Q_vec[k]);
         amrex::Real rho_tot_hi = Rhom_vec[k] * (1. + Q_vec[k]);
-        amrex::Real F = Pm_vec[k] + 0.5*rho_tot_hi*CONST_GRAV*dz + C;
+        amrex::Real F = Pm_vec[k] + 0.5*rho_tot_hi*grav*dz + C;
 
         // Do iterations
         if (std::abs(F)>tol) HSEutils::Newton_Raphson_hse(tol, R_d/Cp_d, dz,
-                                                          CONST_GRAV, C, Thetad_vec[k],
+                                                          grav, C, Thetad_vec[k],
                                                           Q_vec[k], Q_vec[k],
                                                           Pm_vec[k], Rhom_vec[k], F);
     } // k
@@ -178,7 +178,7 @@ calc_rho_p (const int& kmax,
         amrex::Real dz = z_vec[k+1]-z_vec[k];
         Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess.
         for (int it=0; it<maxiter; it++) {
-            Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV;
+            Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*grav;
             if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0;
             Rhod_vec[k] = (p_0/(R_d*Thetam_vec[k]))*std::pow(Pd_vec[k]/p_0, iGamma);
         } // it
diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp
index 9c58ef537..b0f869699 100644
--- a/Source/Initialization/ERF_init_from_metgrid.cpp
+++ b/Source/Initialization/ERF_init_from_metgrid.cpp
@@ -1,7 +1,7 @@
 /**
  * \file ERF_init_from_metgrid.cpp
  */
-
+#include <ERF_Constants.H>
 #include <ERF_Metgrid_utils.H>
 
 using namespace amrex;
@@ -846,6 +846,8 @@ init_base_state_from_metgrid (const bool use_moisture,
     // Make sure this lives on CPU and GPU
     Arena_Used = The_Pinned_Arena();
 #endif
+    // Expose for copy to GPU
+    Real grav = CONST_GRAV;
 
     { // set pressure and density at initialization.
         const Array4<Real>& r_hse_arr  = r_hse_fab.array();
@@ -872,7 +874,7 @@ init_base_state_from_metgrid (const bool use_moisture,
             z_vec[kmax+1] =  new_z(i,j,kmax+1);
 
             calc_rho_p(kmax, flag_psfc_vec[0], orig_psfc(i,j,0),
-                       Thetad_vec, Thetam_vec, Q_vec, z_vec,
+                       grav, Thetad_vec, Thetam_vec, Q_vec, z_vec,
                        Rhod_vec, Rhom_vec, Pd_vec, Pm_vec);
 
             for (int k=0; k<=kmax; k++) {
@@ -976,7 +978,7 @@ init_base_state_from_metgrid (const bool use_moisture,
             z_vec[kmax+1] = new_z(i,j,kmax+1);
 
             calc_rho_p(kmax, flag_psfc_vec[it], orig_psfc(i,j,0),
-                       Thetad_vec, Thetam_vec, Q_vec, z_vec,
+                       grav, Thetad_vec, Thetam_vec, Q_vec, z_vec,
                        Rhod_vec, Rhom_vec, Pd_vec, Pm_vec);
 
             for (int k=0; k<=kmax; k++) {