diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst
index 52938190e..b8f033ab1 100644
--- a/Docs/sphinx_doc/Inputs.rst
+++ b/Docs/sphinx_doc/Inputs.rst
@@ -771,7 +771,7 @@ and "Centered_6th".
 The allowed advection types for the dry and moist scalars are
 "Centered_2nd", "Upwind_3rd", "Blended_3rd4th", "Centered_4th", "Upwind_5th", "Blended_5th6th",
 "Centered_6th" and in addition,
-"WENO3", "WENOZ3", "WENOMZQ3", "WENO5", and "WENOZ5."
+"WENO3", "WENOZ3", "WENOMZQ3", "WENO5", "WENOZ5", "WENO7", and "WENOZ7."
 
 Note: if using WENO schemes, the horizontal and vertical advection types must be set to
 the same string.
diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_WENO b/Exec/RegTests/ScalarAdvDiff/inputs_WENO
index 834a4508a..8bd32f343 100644
--- a/Exec/RegTests/ScalarAdvDiff/inputs_WENO
+++ b/Exec/RegTests/ScalarAdvDiff/inputs_WENO
@@ -19,10 +19,10 @@ erf.cfl            = 0.5     # cfl number for hyperbolic system
 
 erf.dycore_horiz_adv_type    = "Centered_2nd"
 erf.dycore_vert_adv_type     = "Centered_2nd"
-erf.dryscal_horiz_adv_type   = "WENO5"
-erf.dryscal_vert_adv_type    = "WENO5"
-erf.moistscal_horiz_adv_type = "WENO5"
-erf.moistscal_vert_adv_type  = "WENO5"
+erf.dryscal_horiz_adv_type   = "WENO7"
+erf.dryscal_vert_adv_type    = "WENO7"
+erf.moistscal_horiz_adv_type = "WENO7"
+erf.moistscal_vert_adv_type  = "WENO7"
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval   = 1       # timesteps between computing mass
diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z b/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z
index 70d3485ef..774bb4e00 100644
--- a/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z
+++ b/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z
@@ -19,11 +19,11 @@ erf.cfl            = 0.5     # cfl number for hyperbolic system
 
 # TEST WENO SCHEMES
 erf.dycore_horiz_adv_type    = "Centered_2nd"
-erf.dycore_horiz_adv_type    = "Centered_2nd"
-erf.dryscal_vert_adv_type    = "WENO3"
-erf.dryscal_vert_adv_type    = "WENO3"
-erf.moistscal_vert_adv_type  = "WENO3"
-erf.moistscal_vert_adv_type  = "WENO3"
+erf.dycore_vert_adv_type     = "Centered_2nd"
+erf.dryscal_horiz_adv_type   = "WENOZ7"
+erf.dryscal_vert_adv_type    = "WENOZ7"
+erf.moistscal_horiz_adv_type = "WENOZ7"
+erf.moistscal_vert_adv_type  = "WENOZ7"
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval   = 1       # timesteps between computing mass
diff --git a/Source/Advection/ERF_AdvectionSrcForState.cpp b/Source/Advection/ERF_AdvectionSrcForState.cpp
index 74b4f238c..9db3d7c0d 100644
--- a/Source/Advection/ERF_AdvectionSrcForState.cpp
+++ b/Source/Advection/ERF_AdvectionSrcForState.cpp
@@ -235,6 +235,11 @@ AdvectionSrcForScalars (const Real& dt,
                                                        avg_xmom, avg_ymom, avg_zmom,
                                                        horiz_upw_frac, vert_upw_frac);
             break;
+        case AdvType::Weno_7:
+            AdvectionSrcForScalarsWrapper<WENO7,WENO7>(bx, ncomp, icomp, flx_arr, cell_prim,
+                                                       avg_xmom, avg_ymom, avg_zmom,
+                                                       horiz_upw_frac, vert_upw_frac);
+            break;
         case AdvType::Weno_3Z:
             AdvectionSrcForScalarsWrapper<WENO_Z3,WENO_Z3>(bx, ncomp, icomp, flx_arr, cell_prim,
                                                            avg_xmom, avg_ymom, avg_zmom,
@@ -250,6 +255,11 @@ AdvectionSrcForScalars (const Real& dt,
                                                            avg_xmom, avg_ymom, avg_zmom,
                                                            horiz_upw_frac, vert_upw_frac);
             break;
+        case AdvType::Weno_7Z:
+            AdvectionSrcForScalarsWrapper<WENO_Z7,WENO_Z7>(bx, ncomp, icomp, flx_arr, cell_prim,
+                                                           avg_xmom, avg_ymom, avg_zmom,
+                                                           horiz_upw_frac, vert_upw_frac);
+            break;
         default:
             AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
         }
diff --git a/Source/DataStructs/ERF_AdvStruct.H b/Source/DataStructs/ERF_AdvStruct.H
index dfa5bc9f7..5ce564338 100644
--- a/Source/DataStructs/ERF_AdvStruct.H
+++ b/Source/DataStructs/ERF_AdvStruct.H
@@ -120,7 +120,9 @@ struct AdvChoice {
              (dryscal_horiz_adv_string == "WENOZ3"      ) ||
              (dryscal_horiz_adv_string == "WENOMZQ3"    ) ||
              (dryscal_horiz_adv_string == "WENO5"       ) ||
-             (dryscal_horiz_adv_string == "WENOZ5"      ) )
+             (dryscal_horiz_adv_string == "WENOZ5"      ) ||
+             (dryscal_horiz_adv_string == "WENO7"       ) ||
+             (dryscal_horiz_adv_string == "WENOZ7"      ) )
         {
             dryscal_horiz_adv_type = adv_type_convert_string_to_advtype(dryscal_horiz_adv_string);
         }
@@ -134,9 +136,11 @@ struct AdvChoice {
              (dryscal_vert_adv_string == "Centered_6th") ||
              (dryscal_vert_adv_string == "WENO3"       ) ||
              (dryscal_vert_adv_string == "WENOZ3"      ) ||
-             (dryscal_vert_adv_string == "WENOMZQ3"     ) ||
+             (dryscal_vert_adv_string == "WENOMZQ3"    ) ||
              (dryscal_vert_adv_string == "WENO5"       ) ||
-             (dryscal_vert_adv_string == "WENOZ5"      ) )
+             (dryscal_vert_adv_string == "WENOZ5"      ) ||
+             (dryscal_vert_adv_string == "WENO7"       ) ||
+             (dryscal_vert_adv_string == "WENOZ7"      ))
         {
             dryscal_vert_adv_type = adv_type_convert_string_to_advtype(dryscal_vert_adv_string);
         }
@@ -152,7 +156,9 @@ struct AdvChoice {
              (moistscal_horiz_adv_string == "WENOZ3"      ) ||
              (moistscal_horiz_adv_string == "WENOMZQ3"    ) ||
              (moistscal_horiz_adv_string == "WENO5"       ) ||
-             (moistscal_horiz_adv_string == "WENOZ5"      ) )
+             (moistscal_horiz_adv_string == "WENOZ5"      ) ||
+             (moistscal_horiz_adv_string == "WENO7"       ) ||
+             (moistscal_horiz_adv_string == "WENOZ7"      ))
         {
             moistscal_horiz_adv_type = adv_type_convert_string_to_advtype(moistscal_horiz_adv_string);
         }
@@ -168,7 +174,9 @@ struct AdvChoice {
              (moistscal_vert_adv_string == "WENOZ3"      ) ||
              (moistscal_vert_adv_string == "WENOMZQ3"    ) ||
              (moistscal_vert_adv_string == "WENO5"       ) ||
-             (moistscal_vert_adv_string == "WENOZ5"      ) )
+             (moistscal_vert_adv_string == "WENOZ5"      ) ||
+             (moistscal_vert_adv_string == "WENO7"       ) ||
+             (moistscal_vert_adv_string == "WENOZ7"      ))
         {
             moistscal_vert_adv_type = adv_type_convert_string_to_advtype(moistscal_vert_adv_string);
         }
@@ -224,6 +232,10 @@ struct AdvChoice {
             return "WENOZ5";
         } else if (adv_int == AdvType::Weno_3MZQ) {
             return "WENOMZQ3";
+        } else if (adv_int == AdvType::Weno_7) {
+            return "WENO7";
+        } else if (adv_int == AdvType::Weno_7Z) {
+            return "WENOZ7";
         } else {
             return "Unknown";
         }
@@ -251,6 +263,10 @@ struct AdvChoice {
             return AdvType::Weno_5Z;
         } else if (adv_string == "WENOMZQ3") {
             return AdvType::Weno_3MZQ;
+        } else if (adv_string == "WENO7") {
+            return AdvType::Weno_7;
+        } else if (adv_string == "WENOZ7") {
+            return AdvType::Weno_7Z;
         } else {
             return AdvType::Unknown;
         }
diff --git a/Source/ERF.H b/Source/ERF.H
index c69751c69..e2e40b88c 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -1094,6 +1094,16 @@ private:
                || (advChoice.dryscal_horiz_adv_type   == AdvType::Weno_5Z)
                || (advChoice.dryscal_vert_adv_type    == AdvType::Weno_5Z) )
             { return 3; }
+            else if (
+                  (advChoice.dryscal_horiz_adv_type   == AdvType::Weno_7)
+               || (advChoice.dryscal_vert_adv_type    == AdvType::Weno_7)
+               || (advChoice.moistscal_horiz_adv_type == AdvType::Weno_7)
+               || (advChoice.moistscal_vert_adv_type  == AdvType::Weno_7)
+               || (advChoice.moistscal_horiz_adv_type == AdvType::Weno_7Z)
+               || (advChoice.moistscal_vert_adv_type  == AdvType::Weno_7Z)
+               || (advChoice.dryscal_horiz_adv_type   == AdvType::Weno_7Z)
+               || (advChoice.dryscal_vert_adv_type    == AdvType::Weno_7Z) )
+            { return 4; }
             else
             { return 2; }
         }
diff --git a/Source/ERF_IndexDefines.H b/Source/ERF_IndexDefines.H
index 7724e8a36..811946464 100644
--- a/Source/ERF_IndexDefines.H
+++ b/Source/ERF_IndexDefines.H
@@ -193,6 +193,8 @@ enum struct AdvType : int {
     Weno_5        = 108,
     Weno_5Z       = 109,
     Weno_3MZQ     = 110,
-    Unknown       = 111
+    Weno_7        = 111,
+    Weno_7Z       = 112,
+    Unknown       = 113
 };
 #endif
diff --git a/Source/Utils/ERF_Interpolation_WENO.H b/Source/Utils/ERF_Interpolation_WENO.H
index e6325424c..a79fac1ad 100644
--- a/Source/Utils/ERF_Interpolation_WENO.H
+++ b/Source/Utils/ERF_Interpolation_WENO.H
@@ -67,12 +67,12 @@ struct WENO3
     AMREX_FORCE_INLINE
     void
     InterpolateInZ (const int& i,
-                       const int& j,
-                       const int& k,
-                       const int& qty_index,
-                       amrex::Real& val_lo,
-                       amrex::Real upw_lo,
-                       const amrex::Real /*upw_frac*/) const
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
     {
         // Data to interpolate on
         amrex::Real sp1 = m_phi(i  , j  , k+1, qty_index);
@@ -153,7 +153,7 @@ struct WENO5
         if (upw_lo > tol) {
             val_lo = Evaluate(sm3,sm2,sm1,s  ,sp1);
         } else if (upw_lo < -tol) {
-            val_lo = Evaluate(sp2,sp1,s,sm1,sm2);
+            val_lo = Evaluate(sp2,sp1,s  ,sm1,sm2);
         } else {
             val_lo = 0.5 * (s + sm1);
         }
@@ -181,7 +181,7 @@ struct WENO5
         if (upw_lo > tol) {
             val_lo = Evaluate(sm3,sm2,sm1,s  ,sp1);
         } else if (upw_lo < -tol) {
-            val_lo = Evaluate(sp2,sp1,s,sm1,sm2);
+            val_lo = Evaluate(sp2,sp1,s  ,sm1,sm2);
         } else {
             val_lo = 0.5 * (s + sm1);
         }
@@ -191,12 +191,12 @@ struct WENO5
     AMREX_FORCE_INLINE
     void
     InterpolateInZ (const int& i,
-                       const int& j,
-                       const int& k,
-                       const int& qty_index,
-                       amrex::Real& val_lo,
-                       amrex::Real upw_lo,
-                       const amrex::Real /*upw_frac*/) const
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
     {
         // Data to interpolate on
         amrex::Real sp2 = m_phi(i  , j  , k+2, qty_index);
@@ -209,7 +209,7 @@ struct WENO5
         if (upw_lo > tol) {
             val_lo = Evaluate(sm3,sm2,sm1,s  ,sp1);
         } else if (upw_lo < -tol) {
-            val_lo = Evaluate(sp2,sp1,s,sm1,sm2);
+            val_lo = Evaluate(sp2,sp1,s  ,sm1,sm2);
         } else {
             val_lo = 0.5 * (s + sm1);
         }
@@ -258,4 +258,184 @@ private:
     static constexpr amrex::Real g2=(3.0/5.0);
     static constexpr amrex::Real g3=(3.0/10.0);
 };
+
+/**
+ * Interpolation operators used for WENO-7 scheme
+ */
+struct WENO7
+{
+    WENO7 (const amrex::Array4<const amrex::Real>& phi)
+        : m_phi(phi) {}
+
+    AMREX_GPU_DEVICE
+    AMREX_FORCE_INLINE
+    void
+    InterpolateInX (const int& i,
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
+    {
+        // Data to interpolate on
+        amrex::Real sp3 = m_phi(i+3, j  , k  , qty_index);
+        amrex::Real sp2 = m_phi(i+2, j  , k  , qty_index);
+        amrex::Real sp1 = m_phi(i+1, j  , k  , qty_index);
+        amrex::Real s   = m_phi(i  , j  , k  , qty_index);
+        amrex::Real sm1 = m_phi(i-1, j  , k  , qty_index);
+        amrex::Real sm2 = m_phi(i-2, j  , k  , qty_index);
+        amrex::Real sm3 = m_phi(i-3, j  , k  , qty_index);
+        amrex::Real sm4 = m_phi(i-4, j  , k  , qty_index);
+
+        if (upw_lo > tol) {
+            val_lo = Evaluate(sm4,sm3,sm2,sm1,s  ,sp1,sp2);
+        } else if (upw_lo < -tol) {
+            val_lo = Evaluate(sp3,sp2,sp1,s  ,sm1,sm2,sm3);
+        } else {
+            val_lo = 0.5 * (s + sm1);
+        }
+    }
+
+    AMREX_GPU_DEVICE
+    AMREX_FORCE_INLINE
+    void
+    InterpolateInY (const int& i,
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
+    {
+        // Data to interpolate on
+        amrex::Real sp3 = m_phi(i  , j+3, k  , qty_index);
+        amrex::Real sp2 = m_phi(i  , j+2, k  , qty_index);
+        amrex::Real sp1 = m_phi(i  , j+1, k  , qty_index);
+        amrex::Real s   = m_phi(i  , j  , k  , qty_index);
+        amrex::Real sm1 = m_phi(i  , j-1, k  , qty_index);
+        amrex::Real sm2 = m_phi(i  , j-2, k  , qty_index);
+        amrex::Real sm3 = m_phi(i  , j-3, k  , qty_index);
+        amrex::Real sm4 = m_phi(i  , j-4, k  , qty_index);
+
+        if (upw_lo > tol) {
+            val_lo = Evaluate(sm4,sm3,sm2,sm1,s  ,sp1,sp2);
+        } else if (upw_lo < -tol) {
+            val_lo = Evaluate(sp3,sp2,sp1,s  ,sm1,sm2,sm3);
+        } else {
+            val_lo = 0.5 * (s + sm1);
+        }
+    }
+
+    AMREX_GPU_DEVICE
+    AMREX_FORCE_INLINE
+    void
+    InterpolateInZ (const int& i,
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
+    {
+        // Data to interpolate on
+        amrex::Real sp3 = m_phi(i  , j  , k+2, qty_index);
+        amrex::Real sp2 = m_phi(i  , j  , k+2, qty_index);
+        amrex::Real sp1 = m_phi(i  , j  , k+1, qty_index);
+        amrex::Real s   = m_phi(i  , j  , k  , qty_index);
+        amrex::Real sm1 = m_phi(i  , j  , k-1, qty_index);
+        amrex::Real sm2 = m_phi(i  , j  , k-2, qty_index);
+        amrex::Real sm3 = m_phi(i  , j  , k-3, qty_index);
+        amrex::Real sm4 = m_phi(i  , j  , k-3, qty_index);
+
+        if (upw_lo > tol) {
+            val_lo = Evaluate(sm4,sm3,sm2,sm1,s  ,sp1,sp2);
+        } else if (upw_lo < -tol) {
+            val_lo = Evaluate(sp3,sp2,sp1,s  ,sm1,sm2,sm3);
+        } else {
+            val_lo = 0.5 * (s + sm1);
+        }
+    }
+
+    AMREX_GPU_DEVICE
+    AMREX_FORCE_INLINE
+    amrex::Real
+    Evaluate (const amrex::Real& sm3,
+              const amrex::Real& sm2,
+              const amrex::Real& sm1,
+              const amrex::Real& s  ,
+              const amrex::Real& sp1,
+              const amrex::Real& sp2,
+              const amrex::Real& sp3) const
+    {
+        // Smoothing factors
+        amrex::Real b1 = ( sm3 * sm3 * 6649./2880.0
+                         - sm3 * sm2 * 2623./160.0
+                         + sm3 * sm1 * 9449./480.0
+                         - sm3 * s   * 11389./1440.0
+                         + sm2 * sm2 * 28547./960.0
+                         - sm2 * sm1 * 35047./480.0
+                         + sm2 * s   * 14369./480.0
+                         + sm1 * sm1 * 44747./960.0
+                         - sm1 * s   * 6383./160.0
+                         + s   * s   * 25729./2880.0 );
+        amrex::Real b2 = ( sm2 * sm2 * 3169/2880.0
+                         - sm2 * sm1 * 3229/480.0
+                         + sm2 * s   * 3169/480.0
+                         - sm2 * sp1 * 2989/1440.0
+                         + sm1 * sm1 * 11147/960.0
+                         - sm1 * s   * 11767/480.0
+                         + sm1 * sp1 * 1283/160.0
+                         + s   * s   * 13667/960.0
+                         - s   * sp1 * 5069/480.0
+                         + sp1 * sp1 * 6649/2880.0 );
+        amrex::Real b3 = ( sm1 * sm1 * 6649./2880.0
+                         - sm1 * s   * 5069./480.0
+                         + sm1 * sp1 * 1283./160.0
+                         - sm1 * sp2 * 2989./1440.0
+                         + s   * s   * 13667./960.0
+                         - s   * sp1 * 11767./480.0
+                         + s   * sp2 * 3169./480.0
+                         + sp1 * sp1 * 11147./960.0
+                         - sp1 * sp2 * 3229./480.0
+                         + sp2 * sp2 * 3169./2880.0 );
+        amrex::Real b4 = ( s   * s   * 25729./2880.
+                         - s   * sp1 * 6383./160.
+                         + s   * sp2 * 14369./480.
+                         - s   * sp3 * 11389./1440.
+                         + sp1 * sp1 * 44747./960.
+                         - sp1 * sp2 * 35047./480.
+                         + sp1 * sp3 * 9449./480.
+                         + sp2 * sp2 * 28547./960.
+                         - sp2 * sp3 * 2623./160.
+                         + sp3 * sp3 * 6649./2880. );
+
+        // Weight factors
+        amrex::Real w1 = g1 / ( (eps + b1) * (eps + b1) );
+        amrex::Real w2 = g2 / ( (eps + b2) * (eps + b2) );
+        amrex::Real w3 = g3 / ( (eps + b3) * (eps + b3) );
+        amrex::Real w4 = g4 / ( (eps + b4) * (eps + b4) );
+
+        // Weight factor norm
+        amrex::Real wsum = w1 + w2 + w3 + w4;
+
+        // Taylor expansions
+        amrex::Real v1 = (-0.3125)*sm3 + ( 1.3125)*sm2 + (-2.1875)*sm1 + ( 2.1875)*s;
+        amrex::Real v2 = ( 0.0625)*sm2 + (-0.3125)*sm1 + ( 0.9375)*s   + ( 0.3125)*sp1;
+        amrex::Real v3 = (-0.0625)*sm1 + ( 0.5625)*s   + ( 0.5625)*sp1 + (-0.0625)*sp2;
+        amrex::Real v4 = ( 0.3125)*s   + ( 0.9375)*sp1 + (-0.3125)*sp2 + ( 0.0625)*sp3;
+
+        // Interpolated value
+        return ( (w1 * v1 + w2 * v2 + w3 * v3 + w4 * v4) / (wsum) );
+    }
+
+private:
+    amrex::Array4<const amrex::Real> m_phi;   // Quantity to interpolate
+    const amrex::Real eps=1.0e-8;
+    const amrex::Real tol=1.0e-12;
+    static constexpr amrex::Real g1=( 1.0/64.0);
+    static constexpr amrex::Real g2=(21.0/64.0);
+    static constexpr amrex::Real g3=(35.0/64.0);
+    static constexpr amrex::Real g4=( 7.0/64.0);
+};
 #endif
diff --git a/Source/Utils/ERF_Interpolation_WENO_Z.H b/Source/Utils/ERF_Interpolation_WENO_Z.H
index 62e008e1a..bd642d9c6 100644
--- a/Source/Utils/ERF_Interpolation_WENO_Z.H
+++ b/Source/Utils/ERF_Interpolation_WENO_Z.H
@@ -67,12 +67,12 @@ struct WENO_Z3
     AMREX_FORCE_INLINE
     void
     InterpolateInZ (const int& i,
-                       const int& j,
-                       const int& k,
-                       const int& qty_index,
-                       amrex::Real& val_lo,
-                       amrex::Real upw_lo,
-                       const amrex::Real /*upw_frac*/) const
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
     {
         // Data to interpolate on
         amrex::Real sp1 = m_phi(i  , j  , k+1, qty_index);
@@ -188,12 +188,12 @@ struct WENO_MZQ3
     AMREX_FORCE_INLINE
     void
     InterpolateInZ (const int& i,
-                       const int& j,
-                       const int& k,
-                       const int& qty_index,
-                       amrex::Real& val_lo,
-                       amrex::Real upw_lo,
-                       const amrex::Real /*upw_frac*/) const
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
     {
         // Data to interpolate on
         amrex::Real sp1 = m_phi(i  , j  , k+1, qty_index);
@@ -318,12 +318,12 @@ struct WENO_Z5
     AMREX_FORCE_INLINE
     void
     InterpolateInZ (const int& i,
-                       const int& j,
-                       const int& k,
-                       const int& qty_index,
-                       amrex::Real& val_lo,
-                       amrex::Real upw_lo,
-                       const amrex::Real /*upw_frac*/) const
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
     {
         // Data to interpolate on
         amrex::Real sp2 = m_phi(i  , j  , k+2, qty_index);
@@ -386,4 +386,185 @@ private:
     static constexpr amrex::Real g2=(3.0/5.0);
     static constexpr amrex::Real g3=(3.0/10.0);
 };
+
+/**
+ * Interpolation operators used for WENO_Z-7 scheme
+ */
+struct WENO_Z7
+{
+    WENO_Z7 (const amrex::Array4<const amrex::Real>& phi)
+        : m_phi(phi) {}
+
+    AMREX_GPU_DEVICE
+    AMREX_FORCE_INLINE
+    void
+    InterpolateInX (const int& i,
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
+    {
+        // Data to interpolate on
+        amrex::Real sp3 = m_phi(i+3, j  , k  , qty_index);
+        amrex::Real sp2 = m_phi(i+2, j  , k  , qty_index);
+        amrex::Real sp1 = m_phi(i+1, j  , k  , qty_index);
+        amrex::Real s   = m_phi(i  , j  , k  , qty_index);
+        amrex::Real sm1 = m_phi(i-1, j  , k  , qty_index);
+        amrex::Real sm2 = m_phi(i-2, j  , k  , qty_index);
+        amrex::Real sm3 = m_phi(i-3, j  , k  , qty_index);
+        amrex::Real sm4 = m_phi(i-4, j  , k  , qty_index);
+
+        if (upw_lo > tol) {
+            val_lo = Evaluate(sm4,sm3,sm2,sm1,s  ,sp1,sp2);
+        } else if (upw_lo < -tol) {
+            val_lo = Evaluate(sp3,sp2,sp1,s  ,sm1,sm2,sm3);
+        } else {
+            val_lo = 0.5 * (s + sm1);
+        }
+    }
+
+    AMREX_GPU_DEVICE
+    AMREX_FORCE_INLINE
+    void
+    InterpolateInY (const int& i,
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
+    {
+        // Data to interpolate on
+        amrex::Real sp3 = m_phi(i  , j+3, k  , qty_index);
+        amrex::Real sp2 = m_phi(i  , j+2, k  , qty_index);
+        amrex::Real sp1 = m_phi(i  , j+1, k  , qty_index);
+        amrex::Real s   = m_phi(i  , j  , k  , qty_index);
+        amrex::Real sm1 = m_phi(i  , j-1, k  , qty_index);
+        amrex::Real sm2 = m_phi(i  , j-2, k  , qty_index);
+        amrex::Real sm3 = m_phi(i  , j-3, k  , qty_index);
+        amrex::Real sm4 = m_phi(i  , j-4, k  , qty_index);
+
+        if (upw_lo > tol) {
+            val_lo = Evaluate(sm4,sm3,sm2,sm1,s  ,sp1,sp2);
+        } else if (upw_lo < -tol) {
+            val_lo = Evaluate(sp3,sp2,sp1,s  ,sm1,sm2,sm3);
+        } else {
+            val_lo = 0.5 * (s + sm1);
+        }
+    }
+
+    AMREX_GPU_DEVICE
+    AMREX_FORCE_INLINE
+    void
+    InterpolateInZ (const int& i,
+                    const int& j,
+                    const int& k,
+                    const int& qty_index,
+                    amrex::Real& val_lo,
+                    amrex::Real upw_lo,
+                    const amrex::Real /*upw_frac*/) const
+    {
+        // Data to interpolate on
+        amrex::Real sp3 = m_phi(i  , j  , k+2, qty_index);
+        amrex::Real sp2 = m_phi(i  , j  , k+2, qty_index);
+        amrex::Real sp1 = m_phi(i  , j  , k+1, qty_index);
+        amrex::Real s   = m_phi(i  , j  , k  , qty_index);
+        amrex::Real sm1 = m_phi(i  , j  , k-1, qty_index);
+        amrex::Real sm2 = m_phi(i  , j  , k-2, qty_index);
+        amrex::Real sm3 = m_phi(i  , j  , k-3, qty_index);
+        amrex::Real sm4 = m_phi(i  , j  , k-3, qty_index);
+
+        if (upw_lo > tol) {
+            val_lo = Evaluate(sm4,sm3,sm2,sm1,s  ,sp1,sp2);
+        } else if (upw_lo < -tol) {
+            val_lo = Evaluate(sp3,sp2,sp1,s  ,sm1,sm2,sm3);
+        } else {
+            val_lo = 0.5 * (s + sm1);
+        }
+    }
+
+    AMREX_GPU_DEVICE
+    AMREX_FORCE_INLINE
+    amrex::Real
+    Evaluate (const amrex::Real& sm3,
+              const amrex::Real& sm2,
+              const amrex::Real& sm1,
+              const amrex::Real& s  ,
+              const amrex::Real& sp1,
+              const amrex::Real& sp2,
+              const amrex::Real& sp3) const
+    {
+        // Smoothing factors
+        amrex::Real b1 = ( sm3 * sm3 * 6649./2880.0
+                         - sm3 * sm2 * 2623./160.0
+                         + sm3 * sm1 * 9449./480.0
+                         - sm3 * s   * 11389./1440.0
+                         + sm2 * sm2 * 28547./960.0
+                         - sm2 * sm1 * 35047./480.0
+                         + sm2 * s   * 14369./480.0
+                         + sm1 * sm1 * 44747./960.0
+                         - sm1 * s   * 6383./160.0
+                         + s   * s   * 25729./2880.0 );
+        amrex::Real b2 = ( sm2 * sm2 * 3169/2880.0
+                         - sm2 * sm1 * 3229/480.0
+                         + sm2 * s   * 3169/480.0
+                         - sm2 * sp1 * 2989/1440.0
+                         + sm1 * sm1 * 11147/960.0
+                         - sm1 * s   * 11767/480.0
+                         + sm1 * sp1 * 1283/160.0
+                         + s   * s   * 13667/960.0
+                         - s   * sp1 * 5069/480.0
+                         + sp1 * sp1 * 6649/2880.0 );
+        amrex::Real b3 = ( sm1 * sm1 * 6649./2880.0
+                         - sm1 * s   * 5069./480.0
+                         + sm1 * sp1 * 1283./160.0
+                         - sm1 * sp2 * 2989./1440.0
+                         + s   * s   * 13667./960.0
+                         - s   * sp1 * 11767./480.0
+                         + s   * sp2 * 3169./480.0
+                         + sp1 * sp1 * 11147./960.0
+                         - sp1 * sp2 * 3229./480.0
+                         + sp2 * sp2 * 3169./2880.0 );
+        amrex::Real b4 = ( s   * s   * 25729./2880.
+                         - s   * sp1 * 6383./160.
+                         + s   * sp2 * 14369./480.
+                         - s   * sp3 * 11389./1440.
+                         + sp1 * sp1 * 44747./960.
+                         - sp1 * sp2 * 35047./480.
+                         + sp1 * sp3 * 9449./480.
+                         + sp2 * sp2 * 28547./960.
+                         - sp2 * sp3 * 2623./160.
+                         + sp3 * sp3 * 6649./2880. );
+
+        // Weight factors
+        amrex::Real t5 = std::abs(b1 - b2 - b3 + b4);
+        amrex::Real w1 = g1 * ( 1.0 + (t5*t5) / ((eps + b1) * (eps + b1)) );
+        amrex::Real w2 = g2 * ( 1.0 + (t5*t5) / ((eps + b2) * (eps + b2)) );
+        amrex::Real w3 = g3 * ( 1.0 + (t5*t5) / ((eps + b3) * (eps + b3)) );
+        amrex::Real w4 = g4 * ( 1.0 + (t5*t5) / ((eps + b3) * (eps + b3)) );
+
+        // Weight factor norm
+        amrex::Real wsum = w1 + w2 + w3 + w4;
+
+        // Taylor expansions
+        amrex::Real v1 = (-0.3125)*sm3 + ( 1.3125)*sm2 + (-2.1875)*sm1 + ( 2.1875)*s;
+        amrex::Real v2 = ( 0.0625)*sm2 + (-0.3125)*sm1 + ( 0.9375)*s   + ( 0.3125)*sp1;
+        amrex::Real v3 = (-0.0625)*sm1 + ( 0.5625)*s   + ( 0.5625)*sp1 + (-0.0625)*sp2;
+        amrex::Real v4 = ( 0.3125)*s   + ( 0.9375)*sp1 + (-0.3125)*sp2 + ( 0.0625)*sp3;
+
+        // Interpolated value
+        return ( (w1 * v1 + w2 * v2 + w3 * v3 + w4 * v4) / (wsum) );
+    }
+
+private:
+    amrex::Array4<const amrex::Real> m_phi;   // Quantity to interpolate
+    const amrex::Real eps=1.0e-40;
+    const amrex::Real tol=1.0e-12;
+    static constexpr amrex::Real g1=( 1.0/64.0);
+    static constexpr amrex::Real g2=(21.0/64.0);
+    static constexpr amrex::Real g3=(35.0/64.0);
+    static constexpr amrex::Real g4=( 7.0/64.0);
+};
 #endif