diff --git a/CHANGES b/CHANGES
index d12fabc89..76926febd 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,16 @@
+#23.11
+
+  -- AMReX submodule set to 23.11 release hash (ae7b64b)
+
+  -- Bug fix with vertical grid stretching and MOST (#1275)
+
+  -- Bug fix with negative K_turb below surface and MYNN2.5 PBL model (#1271)
+
+  -- Correction to equation of state with moisture (#1263)
+
+  -- Generalized multi-level sponge zones. This allows for nudging of fine solution
+     towards coarse solution with complex box arrays. (#1233)
+
 # 23.10
 
   -- First release since JOSS paper. Too many key ERF PRs to summarize.
diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H
index 1e4204192..c654546f0 100644
--- a/Source/BoundaryConditions/MOSTStress.H
+++ b/Source/BoundaryConditions/MOSTStress.H
@@ -149,7 +149,7 @@ struct adiabatic_charnock
 private:
     most_data mdata;
     similarity_funs sfuns;
-    amrex::Real tol = 1.0e-5;
+    const amrex::Real tol = 1.0e-5;
 };
 
 
@@ -202,7 +202,7 @@ struct adiabatic_mod_charnock
 private:
     most_data mdata;
     similarity_funs sfuns;
-    amrex::Real tol = 1.0e-5;
+    const amrex::Real tol = 1.0e-5;
 };
 
 
@@ -260,7 +260,7 @@ struct surface_flux
 private:
     most_data mdata;
     similarity_funs sfuns;
-    amrex::Real tol = 1.0e-5;
+    const amrex::Real tol = 1.0e-5;
 };
 
 
@@ -323,7 +323,7 @@ struct surface_flux_charnock
 private:
     most_data mdata;
     similarity_funs sfuns;
-    amrex::Real tol = 1.0e-5;
+    const amrex::Real tol = 1.0e-5;
 };
 
 
@@ -387,7 +387,7 @@ struct surface_flux_mod_charnock
 private:
     most_data mdata;
     similarity_funs sfuns;
-    amrex::Real tol = 1.0e-5;
+    const amrex::Real tol = 1.0e-5;
 };
 
 
@@ -447,7 +447,7 @@ struct surface_temp
 private:
     most_data mdata;
     similarity_funs sfuns;
-    amrex::Real tol = 1.0e-5;
+    const amrex::Real tol = 1.0e-5;
 };
 
 
@@ -512,7 +512,7 @@ struct surface_temp_charnock
 private:
     most_data mdata;
     similarity_funs sfuns;
-    amrex::Real tol = 1.0e-5;
+    const amrex::Real tol = 1.0e-5;
 };
 
 
@@ -578,7 +578,7 @@ struct surface_temp_mod_charnock
 private:
     most_data mdata;
     similarity_funs sfuns;
-    amrex::Real tol = 1.0e-5;
+    const amrex::Real tol = 1.0e-5;
 };
 
 
@@ -639,6 +639,7 @@ struct moeng_flux
         rho   = cons_arr(ic,jc,zlo,Rho_comp);
         theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho;
         eta   = eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
+        eta   = amrex::max(eta,eps);
 
         amrex::Real theta_mean  = tm_arr(ic,jc,zlo);
         amrex::Real wsp_mean    = umm_arr(ic,jc,zlo);
@@ -649,7 +650,7 @@ struct moeng_flux
         amrex::Real wsp     = sqrt(velx*velx+vely*vely);
         amrex::Real num1    = (theta-theta_mean)*wsp_mean;
         amrex::Real num2    = (theta_mean-theta_surf)*wsp;
-        amrex::Real moflux  = (std::abs(tstar) > EPS) ?
+        amrex::Real moflux  = (std::abs(tstar) > eps) ?
                               tstar*ustar*(num1+num2)/((theta_mean-theta_surf)*wsp_mean) : 0.0;
         amrex::Real deltaz  = dz * (zlo - k);
 
@@ -700,6 +701,7 @@ struct moeng_flux
                      + cons_arr(ic  ,jc,zlo,Rho_comp) );
         eta   = 0.5 *(  eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
                       + eta_arr(ie  ,je,zlo,EddyDiff::Mom_v) );
+        eta   = amrex::max(eta,eps);
 
         amrex::Real umean    = um_arr(i,j,zlo);
         amrex::Real wsp_mean = 0.5 * ( umm_arr(ic-1,jc,zlo) + umm_arr(ic,jc,zlo) );
@@ -763,6 +765,7 @@ struct moeng_flux
                     + cons_arr(ic,jc  ,zlo,Rho_comp) );
         eta   = 0.5*(  eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
                      + eta_arr(ie,je  ,zlo,EddyDiff::Mom_v) );
+        eta   = amrex::max(eta,eps);
 
         amrex::Real vmean     = vm_arr(i,j,zlo);
         amrex::Real wsp_mean  = 0.5 * ( umm_arr(ic,jc-1,zlo) + umm_arr(ic,jc,zlo) );
@@ -784,7 +787,7 @@ struct moeng_flux
 
 private:
     int zlo;
-    static constexpr amrex::Real EPS = 1e-16;
+    const amrex::Real eps = 1e-16;
 };
 
 
@@ -833,6 +836,7 @@ struct donelan_flux
         rho   = cons_arr(ic,jc,zlo,Rho_comp);
         theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho;
         eta   =  eta_arr(ie,je,zlo,EddyDiff::Theta_v); // == rho * alpha [kg/m^3 * m^2/s]
+        eta   = amrex::max(eta,eps);
 
         amrex::Real Cd = 0.0012;
         amrex::Real wsp_mean    = umm_arr(ic,jc,zlo);
@@ -888,6 +892,7 @@ struct donelan_flux
                      + cons_arr(ic  ,jc,zlo,Rho_comp) );
         eta   = 0.5 *(  eta_arr(ie-1,je,zlo,EddyDiff::Mom_v)
                       + eta_arr(ie  ,je,zlo,EddyDiff::Mom_v) );
+        eta   = amrex::max(eta,eps);
 
         amrex::Real Cd       = 0.001;
         const amrex::Real c  = 7e-5;
@@ -955,6 +960,7 @@ struct donelan_flux
                     + cons_arr(ic,jc  ,zlo,Rho_comp) );
         eta   = 0.5*(  eta_arr(ie,je-1,zlo,EddyDiff::Mom_v)
                      + eta_arr(ie,je  ,zlo,EddyDiff::Mom_v) );
+        eta   = amrex::max(eta,eps);
 
         amrex::Real Cd       = 0.001;
         const amrex::Real c  = 7e-5;
@@ -980,5 +986,6 @@ struct donelan_flux
 
 private:
     int zlo;
+    const amrex::Real eps = 1e-16;
 };
 #endif
diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp
index 4827610ff..bf3709cf4 100644
--- a/Source/Diffusion/ComputeTurbulentViscosity.cpp
+++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp
@@ -193,7 +193,7 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu
     for ( amrex::MFIter mfi(eddyViscosity,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
     {
         Box bxcc   = mfi.tilebox();
-        Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc);
+        Box planex = bxcc; planex.setSmall(0, 1); planex.setBig(0, ngc); planex.grow(1,1);
         Box planey = bxcc; planey.setSmall(1, 1); planey.setBig(1, ngc);
         int i_lo   = bxcc.smallEnd(0); int i_hi = bxcc.bigEnd(0);
         int j_lo   = bxcc.smallEnd(1); int j_hi = bxcc.bigEnd(1);
@@ -202,19 +202,21 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu
 
         const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
 
-        // Extrapolate outside the domain in lateral directions
+        // Extrapolate outside the domain in lateral directions (planex owns corner cells)
         if (i_lo == domain.smallEnd(0)) {
             amrex::ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
-                mu_turb(i_lo-i, j, k, EddyDiff::Mom_h) = mu_turb(i_lo, j, k, EddyDiff::Mom_h);
-                mu_turb(i_lo-i, j, k, EddyDiff::Mom_v) = mu_turb(i_lo, j, k, EddyDiff::Mom_v);
+                int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1));
+                mu_turb(i_lo-i, j, k, EddyDiff::Mom_h) = mu_turb(i_lo, lj, k, EddyDiff::Mom_h);
+                mu_turb(i_lo-i, j, k, EddyDiff::Mom_v) = mu_turb(i_lo, lj, k, EddyDiff::Mom_v);
             });
         }
         if (i_hi == domain.bigEnd(0)) {
             amrex::ParallelFor(planex, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
-                mu_turb(i_hi+i, j, k, EddyDiff::Mom_h) = mu_turb(i_hi, j, k, EddyDiff::Mom_h);
-                mu_turb(i_hi+i, j, k, EddyDiff::Mom_v) = mu_turb(i_hi, j, k, EddyDiff::Mom_v);
+                int lj = amrex::min(amrex::max(j, domain.smallEnd(1)), domain.bigEnd(1));
+                mu_turb(i_hi+i, j, k, EddyDiff::Mom_h) = mu_turb(i_hi, lj, k, EddyDiff::Mom_h);
+                mu_turb(i_hi+i, j, k, EddyDiff::Mom_v) = mu_turb(i_hi, lj, k, EddyDiff::Mom_v);
             });
         }
         if (j_lo == domain.smallEnd(1)) {
@@ -336,9 +338,9 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu
                    mu_turb(i, j, k_hi+k, indx_v) = mu_turb(i, j, k_hi, indx_v);
                  });
                  break;
-           }
-       }
-   }
+            }
+        }
+    }
 }
 
 /**
diff --git a/Submodules/AMReX b/Submodules/AMReX
index 388738dc7..ae7b64bcf 160000
--- a/Submodules/AMReX
+++ b/Submodules/AMReX
@@ -1 +1 @@
-Subproject commit 388738dc788cc4a72c33a5d5d9940ef3a37a76b7
+Subproject commit ae7b64bcf6a4ed36dd03e17357cae83b0a394912