diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 569742c872..1845903780 100755
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -123,6 +123,18 @@ steps:
         command: "julia --color=yes --project=examples examples/hybrid/driver.jl --rad allskywithclear --idealized_h2o true --idealized_clouds true --hyperdiff false --config column --z_max 70e3 --z_elem 70 --dz_bottom 100 --dz_top 10000 --t_end 654days --dt 3hours --dt_save_to_sol 30hours --job_id sphere_single_column_radiative_equilibrium_allsky_idealized_clouds"
         artifact_paths: "sphere_single_column_radiative_equilibrium_allsky_idealized_clouds/*"
 
+      - label: ":computer: single column non-orographic gravity wave parameterization"
+        command: "julia --color=yes --project=examples examples/hybrid/driver.jl --hyperdiff false --config column --non_orographic_gravity_wave true --t_end 1500secs --dt 500secs --job_id single_column_nonorographic_gravity_wave --regression_test true"
+        artifact_paths: "single_column_nonorographic_gravity_wave/*"
+
+      - label: ":computer: non-orographic gravity wave parameterization unit test 3d"
+        command: "julia --color=yes --project=examples examples/hybrid/gravitywave_parameterization/gw_test_3d.jl"
+        artifact_paths: "nonorographic_gravity_wave_test_3d/*"
+
+      - label: ":computer: non-orographic gravity wave parameterization unit test single column"
+        command: "julia --color=yes --project=examples examples/hybrid/gravitywave_parameterization/gw_test_single_column.jl"
+        artifact_paths: "nonorographic_gravity_wave_test_single_column/*"
+
   - group: "MPI Examples"
     steps:
 
diff --git a/examples/hybrid/cli_options.jl b/examples/hybrid/cli_options.jl
index 3b06577e63..dcbde45c22 100644
--- a/examples/hybrid/cli_options.jl
+++ b/examples/hybrid/cli_options.jl
@@ -226,6 +226,10 @@ function parse_commandline()
         help = "Save most of the tc aux state to HDF5 file [`false` (default), `true`]"
         arg_type = Bool
         default = false
+        "--non_orographic_gravity_wave"
+        help = "Apply parameterization for convective gravity wave forcing on horizontal mean flow"
+        arg_type = Bool
+        default = false
     end
     parsed_args = ArgParse.parse_args(ARGS, s)
     return (s, parsed_args)
diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl
index 884217723a..8dc848ce4c 100644
--- a/examples/hybrid/driver.jl
+++ b/examples/hybrid/driver.jl
@@ -157,6 +157,8 @@ function additional_cache(Y, params, model_spec, dt; use_tempest_mode = false)
             diffuse_momentum,
             coupled,
         ) : NamedTuple(),
+        model_spec.non_orographic_gravity_wave ? gravity_wave_cache(Y, FT) :
+        NamedTuple(),
         (;
             tendency_knobs = (;
                 hs_forcing = forcing_type isa HeldSuarezForcing,
@@ -166,6 +168,7 @@ function additional_cache(Y, params, model_spec, dt; use_tempest_mode = false)
                 rayleigh_sponge,
                 viscous_sponge,
                 hyperdiff,
+                non_orographic_gravity_wave = model_spec.non_orographic_gravity_wave,
                 has_turbconv = !isnothing(turbconv_model),
             )
         ),
@@ -182,6 +185,7 @@ function additional_tendency!(Yₜ, Y, p, t)
     (; rad_flux, vert_diff, hs_forcing) = p.tendency_knobs
     (; microphy_0M, hyperdiff, has_turbconv) = p.tendency_knobs
     (; rayleigh_sponge, viscous_sponge) = p.tendency_knobs
+    (; non_orographic_gravity_wave) = p.tendency_knobs
     hyperdiff && hyperdiffusion_tendency!(Yₜ, Y, p, t)
     rayleigh_sponge && rayleigh_sponge_tendency!(Yₜ, Y, p, t)
     viscous_sponge && viscous_sponge_tendency!(Yₜ, Y, p, t)
@@ -190,6 +194,7 @@ function additional_tendency!(Yₜ, Y, p, t)
     microphy_0M && zero_moment_microphysics_tendency!(Yₜ, Y, p, t)
     rad_flux && rrtmgp_model_tendency!(Yₜ, Y, p, t)
     has_turbconv && TCU.sgs_flux_tendency!(Yₜ, Y, p, t)
+    non_orographic_gravity_wave && gravity_wave_tendency!(Yₜ, Y, p, t)
 end
 
 ################################################################################
@@ -204,6 +209,12 @@ if parsed_args["trunc_stack_traces"]
 end
 
 include(joinpath("sphere", "baroclinic_wave_utilities.jl"))
+include(
+    joinpath(
+        "gravitywave_parameterization",
+        "gravity_wave_parameterization.jl",
+    ),
+)
 
 
 import ClimaCore: enable_threading
diff --git a/examples/hybrid/gravitywave_parameterization/gravity_wave_parameterization.jl b/examples/hybrid/gravitywave_parameterization/gravity_wave_parameterization.jl
new file mode 100644
index 0000000000..ec7fc61e99
--- /dev/null
+++ b/examples/hybrid/gravitywave_parameterization/gravity_wave_parameterization.jl
@@ -0,0 +1,240 @@
+function gravity_wave_cache(
+    Y,
+    ::Type{FT};
+    source_height = FT(15000),
+    Bm = FT(1.2),
+    F_S0 = FT(4e-3),
+    dc = FT(0.6),
+    cmax = FT(99.6),
+    c0 = FT(0),
+    kwv = FT(2π / 100e5),
+    cw = FT(40.0),
+) where {FT}
+
+    nc = Int(floor(FT(2 * cmax / dc + 1)))
+    c = [FT((n - 1) * dc - cmax) for n in 1:nc]
+
+    return (;
+        gw_source_height = source_height,
+        gw_F_S0 = F_S0,
+        gw_Bm = Bm,
+        gw_c = c,
+        gw_cw = cw,
+        gw_c0 = c0,
+        gw_nk = length(kwv),
+        gw_k = kwv,
+        gw_k2 = kwv .^ 2,
+        ᶜbuoyancy_frequency = similar(Y.c.ρ),
+        ᶜdTdz = similar(Y.c.ρ),
+    )
+end
+
+function gravity_wave_tendency!(Yₜ, Y, p, t)
+    #unpack
+    (;
+        gw_source_height,
+        gw_Bm,
+        gw_F_S0,
+        gw_c,
+        gw_cw,
+        gw_c0,
+        gw_nk,
+        gw_k,
+        gw_k2,
+    ) = p
+    (; ᶜts, ᶜT, ᶜdTdz, ᶜbuoyancy_frequency, params) = p
+    ᶜρ = Y.c.ρ
+    # parameters
+    thermo_params = CAP.thermodynamics_params(params)
+    grav = FT(CAP.grav(params))
+
+    # compute buoyancy frequency
+    @. ᶜT = TD.air_temperature(thermo_params, ᶜts)
+
+    parent(ᶜdTdz) .= parent(Geometry.WVector.(ᶜgradᵥ.(ᶠinterp.(ᶜT))))
+
+    ᶜbuoyancy_frequency =
+        @. (grav / ᶜT) * (ᶜdTdz + grav / TD.cp_m(thermo_params, ᶜts))
+    ᶜbuoyancy_frequency = @. ifelse(
+        ᶜbuoyancy_frequency < FT(2.5e-5),
+        FT(sqrt(2.5e-5)),
+        sqrt(abs(ᶜbuoyancy_frequency)),
+    ) # to avoid small numbers
+    # alternative
+    # ᶜbuoyancy_frequency = [i > 2.5e-5 ? sqrt(i) : sqrt(2.5e-5)  for i in ᶜbuoyancy_frequency]
+    # TODO: create an extra layer at model top so that the gravity wave forcing
+    # .     occurring between the topmost model level and the upper boundary
+    # .     may be calculated
+    #
+
+    # source level: get the index of the level that is closest to the source height (GFDL uses the fist level below instead)
+    source_level = argmin(
+        abs.(
+            unique(parent(Fields.coordinate_field(Y.c).z) .- gw_source_height)
+        ),
+    )
+    ᶠz = Fields.coordinate_field(Y.f).z
+
+    # TODO: prepare input variables for gravity_wave_forcing()
+    #       1. grab column wise data
+    #       2. convert coviant uv to physical uv
+    # need to make everything an array befere
+    u_phy = Geometry.UVVector.(Y.c.uₕ).components.data.:1
+    v_phy = Geometry.UVVector.(Y.c.uₕ).components.data.:2
+
+    uforcing = ones(axes(u_phy))
+    vforcing = similar(v_phy)
+    # bycolume
+    Fields.bycolumn(axes(ᶜρ)) do colidx
+        parent(uforcing[colidx]) .= gravity_wave_forcing(
+            parent(u_phy[colidx]),
+            source_level,
+            gw_F_S0,
+            gw_Bm,
+            gw_c,
+            gw_cw,
+            gw_c0,
+            gw_nk,
+            gw_k,
+            gw_k2,
+            parent(ᶜbuoyancy_frequency[colidx]),
+            parent(ᶜρ[colidx]),
+            parent(ᶠz[colidx]),
+        )
+        parent(vforcing[colidx]) .= gravity_wave_forcing(
+            parent(v_phy[colidx]),
+            source_level,
+            gw_F_S0,
+            gw_Bm,
+            gw_c,
+            gw_cw,
+            gw_c0,
+            gw_nk,
+            gw_k,
+            gw_k2,
+            parent(ᶜbuoyancy_frequency[colidx]),
+            parent(ᶜρ[colidx]),
+            parent(ᶠz[colidx]),
+        )
+
+    end
+
+    @. Yₜ.c.uₕ +=
+        Geometry.Covariant12Vector.(Geometry.UVVector.(uforcing, vforcing))
+
+end
+
+function gravity_wave_forcing(
+    ᶜu,
+    source_level,
+    F_S0,
+    Bm,
+    c,
+    cw,
+    c0,
+    nk,
+    kwv,
+    k2,
+    ᶜbf,
+    ᶜρ,
+    ᶠz,
+)
+
+    nc = length(c)
+
+    # define wave momentum flux (B0) at source level for each phase
+    # speed n, and the sum over all phase speeds (Bsum), which is needed
+    # to calculate the intermittency.
+
+    c_hat0 = c .- ᶜu[source_level]
+    # In GFDL code, flag is always 1: c = c0 * flag + c_hat0 * (1 - flag)
+    Bexp = @. exp(-log(2.0) * ((c - c0) / cw)^2)
+    B0 = @. sign(c_hat0) * Bm * Bexp
+    # In GFDL code, it is: B0 = @. sign(c_hat0) * (Bw * Bexp + Bn * Bexp)
+    # where Bw = Bm is the wide band and Bn is the narrow band
+    Bsum = sum(abs.(B0))
+    if (Bsum == 0.0)
+        error("zero flux input at source level")
+    end
+    # define the intermittency factor eps. Assumes constant Δc.
+    eps = (F_S0 / ᶜρ[source_level] / nk) / Bsum # in GFDL code: eps = (ampl * 1.5 / nk) / Bsum, ampl is equivalent to F_S)/ρ_source 
+
+    ᶜdz = ᶠz[2:end] - ᶠz[1:(end - 1)]
+    wave_forcing = zeros(nc)
+    gwf = zeros(length(ᶜu))
+    for ink in 1:nk # loop over all wave lengths
+
+        mask = ones(nc)  # mask to determine which waves propagate upward
+        for k in source_level:length(ᶜu)
+            fac = FT(0.5) * (ᶜρ[k] / ᶜρ[source_level]) * kwv[ink] / ᶜbf[k]
+
+            ᶜHb = ᶜdz[k] / log(ᶜρ[k - 1] / ᶜρ[k])  # density scale height
+            alp2 = 0.25 / (ᶜHb * ᶜHb)
+            ω_r = sqrt((ᶜbf[k] * ᶜbf[k] * k2[ink]) / (k2[ink] + alp2)) # ω_r (critical frequency that marks total internal reflection)
+
+            fm = FT(0)
+            for n in 1:nc
+                # check only those waves which are still propagating, i.e., mask = 1.0
+                if (mask[n]) == 1.0
+                    c_hat = c[n] - ᶜu[k]
+                    # f phase speed matches the wind speed, remove c(n) from the set of propagating waves.
+                    if c_hat == 0.0
+                        mask[n] = 0.0
+                    else
+                        # define the criterion which determines if wave is reflected at this level (test).
+                        test = abs(c_hat) * kwv[ink] - ω_r
+                        if test >= 0.0
+                            # wave has undergone total internal reflection. remove it from the propagating set.
+                            mask[n] = 0.0
+                        else
+                            # if wave is not reflected at this level, determine if it is
+                            # breaking at this level (Foc >= 0), or if wave speed relative to
+                            # windspeed has changed sign from its value at the source level
+                            # (c_hat0[n] * c_hat <= 0). if it is above the source level and is
+                            # breaking, then add its momentum flux to the accumulated sum at
+                            # this level.
+                            # set mask=0.0 to remove phase speed band c[n] from the set of active 
+                            # waves moving upwards to the next level.
+                            if c_hat0[n] * c_hat <= 0.0
+                                mask[n] = 0.0
+                                if k > source_level
+                                    fm = fm + B0[n]
+                                end
+                            else
+                                Foc = B0[n] / (c_hat)^3 - fac
+                                if Foc >= 0.0
+                                    mask[n] = 0.0
+                                    if k > source_level
+                                        fm = fm + B0[n]
+                                    end
+                                end
+                            end
+                        end # (test >= 0.0)
+                    end #(c_hat == 0.0)
+                end # mask = 0
+
+            end # phase speed loop
+
+            # TODO: GFDL option to dump remaining flux at the top of the model
+
+            # compute the gravity wave momentum flux forcing 
+            # obtained across the entire wave spectrum at this level.
+            if k > source_level
+                rbh = sqrt(ᶜρ[k] * ᶜρ[k - 1])
+                wave_forcing[k] = (ᶜρ[source_level] / rbh) * fm * eps / ᶜdz[k]
+                wave_forcing[k - 1] =
+                    0.5 * (wave_forcing[k - 1] + wave_forcing[k])
+            else
+                wave_forcing[k] = 0.0
+            end
+
+        end # k
+
+        for k in source_level:length(ᶜu)
+            gwf[k] = gwf[k] + wave_forcing[k]
+        end
+
+    end # ink 
+
+    return gwf
+end
diff --git a/examples/hybrid/gravitywave_parameterization/gw_test_3d.jl b/examples/hybrid/gravitywave_parameterization/gw_test_3d.jl
new file mode 100644
index 0000000000..0e6a430b5e
--- /dev/null
+++ b/examples/hybrid/gravitywave_parameterization/gw_test_3d.jl
@@ -0,0 +1,170 @@
+using NCDatasets
+using Dates
+using Interpolations
+using Statistics
+using Plots
+include("gravity_wave_parameterization.jl")
+const FT = Float64
+
+# test Figure 8 of the Alexander and Dunkerton (1999) paper:
+# https://journals.ametsoc.org/view/journals/atsc/56/24/1520-0469_1999_056_4167_aspomf_2.0.co_2.xml?tab_body=pdf
+
+face_z = 0:1e3:0.47e5
+center_z = 0.5 .* (face_z[1:(end - 1)] .+ face_z[2:end])
+
+# compute the source parameters
+function gravity_wave_cache(
+    ::Type{FT};
+    source_height = FT(15000),
+    Bm = FT(1.2),
+    F_S0 = FT(4e-3),
+    dc = FT(0.6),
+    cmax = FT(99.6),
+    c0 = FT(0),
+    kwv = FT(2π / 100e5),
+    cw = FT(40.0),
+) where {FT}
+
+    nc = Int(floor(FT(2 * cmax / dc + 1)))
+    c = [FT((n - 1) * dc - cmax) for n in 1:nc]
+
+    return (;
+        gw_source_height = source_height,
+        gw_F_S0 = F_S0,
+        gw_Bm = Bm,
+        gw_c = c,
+        gw_cw = cw,
+        gw_c0 = c0,
+        gw_nk = length(kwv),
+        gw_k = kwv,
+        gw_k2 = kwv .^ 2,
+    )
+end
+
+params = gravity_wave_cache(FT; Bm = 0.4, cmax = 150, kwv = 2π / 100e3)
+source_level = argmin(abs.(center_z .- params.gw_source_height))
+
+# ERA5 data 1973 Jan: downloaded from Caltech box
+using Pkg.Artifacts
+using ArtifactWrappers
+
+function era_dataset_path()
+    era_dataset = ArtifactWrapper(
+        @__DIR__,
+        isempty(get(ENV, "CI", "")),
+        "era-global",
+        ArtifactFile[ArtifactFile(
+            url = "https://caltech.box.com/shared/static/2489yvlwhnxsrl05jvqnch7fcrd0k1vw.nc",
+            filename = "box-era5-monthly.nc",
+        ),],
+    )
+    return get_data_folder(era_dataset)
+end
+era_data = joinpath(era_dataset_path(), "box-era5-monthly.nc")
+
+nt = NCDataset(era_data) do ds
+    lon = ds["longitude"][:]
+    lat = ds["latitude"][:]
+    lev = ds["level"][:] .* 100
+    time = ds["time"][:]
+    gZ = ds["z"][:]
+    T = ds["t"][:]
+    u = ds["u"][:]
+    (; lon, lat, lev, time, gZ, T, u)
+end
+(; lon, lat, lev, time, gZ, T, u) = nt
+
+# compute density and buoyancy frequency
+R_d = 287.0
+grav = 9.8
+cp_d = 1004.0
+
+Z = gZ ./ grav
+ρ = ones(size(T)) .* reshape(lev, (1, 1, length(lev), 1)) ./ T / R_d
+
+dTdz = zeros(size(T))
+@. dTdz[:, :, 1, :] =
+    (T[:, :, 2, :] - T[:, :, 1, :]) / (Z[:, :, 2, :] - Z[:, :, 1, :])
+@. dTdz[:, :, end, :] =
+    (T[:, :, end, :] - T[:, :, end - 1, :]) /
+    (Z[:, :, end, :] - Z[:, :, end - 1, :])
+@. dTdz[:, :, 2:(end - 1), :] =
+    (T[:, :, 3:end, :] - T[:, :, 1:(end - 2), :]) /
+    (Z[:, :, 3:end, :] - Z[:, :, 1:(end - 2), :])
+bf = @. (grav / T) * (dTdz + grav / cp_d)
+bf = @. ifelse(bf < 2.5e-5, sqrt(2.5e-5), sqrt(abs(bf)))
+
+# interpolation to center_z grid
+center_u = zeros(length(lon), length(lat), length(center_z), length(time))
+center_bf = zeros(length(lon), length(lat), length(center_z), length(time))
+center_ρ = zeros(length(lon), length(lat), length(center_z), length(time))
+for i in 1:length(lon)
+    for j in 1:length(lat)
+        for it in 1:length(time)
+            interp_linear = LinearInterpolation(
+                Z[i, j, :, it][end:-1:1],
+                u[i, j, :, it][end:-1:1],
+                extrapolation_bc = Line(),
+            )
+            center_u[i, j, :, it] = interp_linear.(center_z)
+
+            interp_linear = LinearInterpolation(
+                Z[i, j, :, it][end:-1:1],
+                bf[i, j, :, it][end:-1:1],
+                extrapolation_bc = Line(),
+            )
+            center_bf[i, j, :, it] = interp_linear.(center_z)
+
+            interp_linear = LinearInterpolation(
+                Z[i, j, :, it][end:-1:1],
+                ρ[i, j, :, it][end:-1:1],
+                extrapolation_bc = Line(),
+            )
+            center_ρ[i, j, :, it] = interp_linear.(center_z)
+        end
+    end
+end
+
+# compute zonal mean profile first and apply parameterization 
+center_u_zonalave = mean(center_u, dims = 1)[1, :, :, :]
+center_bf_zonalave = mean(center_bf, dims = 1)[1, :, :, :]
+center_ρ_zonalave = mean(center_ρ, dims = 1)[1, :, :, :]
+
+# Jan
+month = Dates.month.(time)
+
+Jan_u = mean(center_u_zonalave[:, :, month .== 1], dims = 3)[:, :, 1]
+Jan_bf = mean(center_bf_zonalave[:, :, month .== 1], dims = 3)[:, :, 1]
+Jan_ρ = mean(center_ρ_zonalave[:, :, month .== 1], dims = 3)[:, :, 1]
+Jan_uforcing = zeros(length(lat), length(center_z))
+for j in 1:length(lat)
+    Jan_uforcing[j, :] = gravity_wave_forcing(
+        Jan_u[j, :],
+        source_level,
+        params.gw_F_S0,
+        params.gw_Bm,
+        params.gw_c,
+        params.gw_cw,
+        params.gw_c0,
+        params.gw_nk,
+        params.gw_k,
+        params.gw_k2,
+        Jan_bf[j, :],
+        Jan_ρ[j, :],
+        face_z,
+    )
+end
+
+ENV["GKSwstype"] = "nul"
+output_dir = "nonorographic_gravity_wave_test_3d"
+mkpath(output_dir)
+png(
+    contourf(
+        lat[end:-1:1],
+        center_z[source_level:end],
+        86400 * Jan_uforcing[end:-1:1, source_level:end]',
+        color = :balance,
+        clim = (-1, 1),
+    ),
+    joinpath(output_dir, "test-fig8.png"),
+)
diff --git a/examples/hybrid/gravitywave_parameterization/gw_test_single_column.jl b/examples/hybrid/gravitywave_parameterization/gw_test_single_column.jl
new file mode 100644
index 0000000000..db1e97f791
--- /dev/null
+++ b/examples/hybrid/gravitywave_parameterization/gw_test_single_column.jl
@@ -0,0 +1,233 @@
+using NCDatasets
+using Dates
+using Interpolations
+using Statistics
+using Plots
+include("gravity_wave_parameterization.jl")
+
+const FT = Float64
+# single column test Figure 6 of the Alexander and Dunkerton (1999) paper:
+# https://journals.ametsoc.org/view/journals/atsc/56/24/1520-0469_1999_056_4167_aspomf_2.0.co_2.xml?tab_body=pdf
+# zonal mean monthly wind, temperature 1958-1973; at 40N in latitude for Jan, April, July, Oct.
+
+face_z = FT.(0:1e3:0.5e5)
+center_z = FT(0.5) .* (face_z[1:(end - 1)] .+ face_z[2:end])
+
+# compute the source parameters
+function gravity_wave_cache(
+    ::Type{FT};
+    source_height = FT(15000),
+    Bm = FT(1.2),
+    F_S0 = FT(4e-3),
+    dc = FT(0.6),
+    cmax = FT(99.6),
+    c0 = FT(0),
+    kwv = FT(2π / 100e5),
+    cw = FT(40.0),
+) where {FT}
+
+    nc = Int(floor(FT(2 * cmax / dc + 1)))
+    c = [FT((n - 1) * dc - cmax) for n in 1:nc]
+
+    return (;
+        gw_source_height = source_height,
+        gw_F_S0 = F_S0,
+        gw_Bm = Bm,
+        gw_c = c,
+        gw_cw = cw,
+        gw_c0 = c0,
+        gw_nk = length(kwv),
+        gw_k = kwv,
+        gw_k2 = kwv .^ 2,
+    )
+end
+
+params = gravity_wave_cache(FT; Bm = 0.4, cmax = 150, kwv = 2π / 100e3)
+source_level = argmin(abs.(center_z .- params.gw_source_height))
+
+# download ERA5 nc data: wind, temperature, geopotential height
+using Pkg.Artifacts
+using ArtifactWrappers
+
+function era_dataset_path()
+    era_dataset = ArtifactWrapper(
+        @__DIR__,
+        isempty(get(ENV, "CI", "")),
+        "era-single-column",
+        ArtifactFile[ArtifactFile(
+            url = "https://caltech.box.com/shared/static/of5wi39o643a333yy9vbx5pnf0za503g.nc",
+            filename = "box-single_column_test.nc",
+        ),],
+    )
+    return get_data_folder(era_dataset)
+end
+era_data = joinpath(era_dataset_path(), "box-single_column_test.nc")
+
+nt = NCDataset(era_data) do ds
+    # Dimensions:  longitude × latitude × level × time
+    # data at 40N, all longitude, all levels, all time
+    lon = ds["longitude"][:]
+    lat = ds["latitude"][:]
+    lev = ds["level"][:] .* 100
+    time = ds["time"][:]
+    gZ = ds["z"][:, 5, :, :]
+    T = ds["t"][:, 5, :, :]
+    u = ds["u"][:, 5, :, :]
+    (; lon, lat, lev, time, gZ, T, u)
+end
+(; lon, lat, lev, time, gZ, T, u) = nt
+
+# compute density and buoyancy frequency
+R_d = 287.0
+grav = 9.8
+cp_d = 1004.0
+
+Z = gZ ./ grav
+ρ = ones(size(T)) .* reshape(lev, (1, length(lev), 1)) ./ T / R_d
+
+dTdz = zeros(size(T))
+@. dTdz[:, 1, :] = (T[:, 2, :] - T[:, 1, :]) / (Z[:, 2, :] - Z[:, 1, :])
+@. dTdz[:, end, :] =
+    (T[:, end, :] - T[:, end - 1, :]) / (Z[:, end, :] - Z[:, end - 1, :])
+@. dTdz[:, 2:(end - 1), :] =
+    (T[:, 3:end, :] - T[:, 1:(end - 2), :]) /
+    (Z[:, 3:end, :] - Z[:, 1:(end - 2), :])
+bf = @. (grav / T) * (dTdz + grav / cp_d)
+bf = @. ifelse(bf < 2.5e-5, sqrt(2.5e-5), sqrt(abs(bf)))
+
+# interpolation to center_z grid
+center_u = zeros(length(lon), length(center_z), length(time))
+center_bf = zeros(length(lon), length(center_z), length(time))
+center_ρ = zeros(length(lon), length(center_z), length(time))
+for i in 1:length(lon)
+    for it in 1:length(time)
+        interp_linear = LinearInterpolation(
+            Z[i, :, it][end:-1:1],
+            u[i, :, it][end:-1:1],
+            extrapolation_bc = Line(),
+        )
+        center_u[i, :, it] = interp_linear.(center_z)
+
+        interp_linear = LinearInterpolation(
+            Z[i, :, it][end:-1:1],
+            bf[i, :, it][end:-1:1],
+            extrapolation_bc = Line(),
+        )
+        center_bf[i, :, it] = interp_linear.(center_z)
+
+        interp_linear = LinearInterpolation(
+            Z[i, :, it][end:-1:1],
+            ρ[i, :, it][end:-1:1],
+            extrapolation_bc = Line(),
+        )
+        center_ρ[i, :, it] = interp_linear.(center_z)
+    end
+end
+
+# zonal mean
+center_u_mean = mean(center_u, dims = 1)[1, :, :]
+center_bf_mean = mean(center_bf, dims = 1)[1, :, :]
+center_ρ_mean = mean(center_ρ, dims = 1)[1, :, :]
+
+# monthly ave Jan, April, July, Oct
+month = Dates.month.(time)
+
+ENV["GKSwstype"] = "nul"
+output_dir = "nonorographic_gravity_wave_test_single_column"
+mkpath(output_dir)
+
+# Jan
+Jan_u = mean(center_u_mean[:, month .== 1], dims = 2)[:, 1]
+Jan_bf = mean(center_bf_mean[:, month .== 1], dims = 2)[:, 1]
+Jan_ρ = mean(center_ρ_mean[:, month .== 1], dims = 2)[:, 1]
+Jan_uforcing = gravity_wave_forcing(
+    Jan_u,
+    source_level,
+    params.gw_F_S0,
+    params.gw_Bm,
+    params.gw_c,
+    params.gw_cw,
+    params.gw_c0,
+    params.gw_nk,
+    params.gw_k,
+    params.gw_k2,
+    Jan_bf,
+    Jan_ρ,
+    face_z,
+)
+png(
+    plot(Jan_uforcing[source_level:end] * 86400, center_z[source_level:end]),
+    joinpath(output_dir, "fig6jan.png"),
+)
+
+# April
+April_u = mean(center_u_mean[:, month .== 4], dims = 2)[:, 1]
+April_bf = mean(center_bf_mean[:, month .== 4], dims = 2)[:, 1]
+April_ρ = mean(center_ρ_mean[:, month .== 4], dims = 2)[:, 1]
+April_uforcing = gravity_wave_forcing(
+    April_u,
+    source_level,
+    params.gw_F_S0,
+    params.gw_Bm,
+    params.gw_c,
+    params.gw_cw,
+    params.gw_c0,
+    params.gw_nk,
+    params.gw_k,
+    params.gw_k2,
+    April_bf,
+    April_ρ,
+    face_z,
+)
+png(
+    plot(April_uforcing[source_level:end] * 86400, center_z[source_level:end]),
+    joinpath(output_dir, "fig6apr.png"),
+)
+
+# July
+July_u = mean(center_u_mean[:, month .== 7], dims = 2)[:, 1]
+July_bf = mean(center_bf_mean[:, month .== 7], dims = 2)[:, 1]
+July_ρ = mean(center_ρ_mean[:, month .== 7], dims = 2)[:, 1]
+July_uforcing = gravity_wave_forcing(
+    July_u,
+    source_level,
+    params.gw_F_S0,
+    params.gw_Bm,
+    params.gw_c,
+    params.gw_cw,
+    params.gw_c0,
+    params.gw_nk,
+    params.gw_k,
+    params.gw_k2,
+    July_bf,
+    July_ρ,
+    face_z,
+)
+png(
+    plot(July_uforcing[source_level:end] * 86400, center_z[source_level:end]),
+    joinpath(output_dir, "fig6jul.png"),
+)
+
+# Oct
+Oct_u = mean(center_u_mean[:, month .== 10], dims = 2)[:, 1]
+Oct_bf = mean(center_bf_mean[:, month .== 10], dims = 2)[:, 1]
+Oct_ρ = mean(center_ρ_mean[:, month .== 10], dims = 2)[:, 1]
+Oct_uforcing = gravity_wave_forcing(
+    Oct_u,
+    source_level,
+    params.gw_F_S0,
+    params.gw_Bm,
+    params.gw_c,
+    params.gw_cw,
+    params.gw_c0,
+    params.gw_nk,
+    params.gw_k,
+    params.gw_k2,
+    Oct_bf,
+    Oct_ρ,
+    face_z,
+)
+png(
+    plot(Oct_uforcing[source_level:end] * 86400, center_z[source_level:end]),
+    joinpath(output_dir, "fig6oct.png"),
+)
diff --git a/examples/hybrid/staggered_nonhydrostatic_model.jl b/examples/hybrid/staggered_nonhydrostatic_model.jl
index 17a3484066..7f64b65e69 100644
--- a/examples/hybrid/staggered_nonhydrostatic_model.jl
+++ b/examples/hybrid/staggered_nonhydrostatic_model.jl
@@ -50,6 +50,7 @@ const ᶠgradᵥ = Operators.GradientC2F(
     bottom = Operators.SetGradient(Geometry.Covariant3Vector(FT(0))),
     top = Operators.SetGradient(Geometry.Covariant3Vector(FT(0))),
 )
+const ᶜgradᵥ = Operators.GradientF2C()
 const ᶠcurlᵥ = Operators.CurlC2F(
     bottom = Operators.SetCurl(Geometry.Contravariant12Vector(FT(0), FT(0))),
     top = Operators.SetCurl(Geometry.Contravariant12Vector(FT(0), FT(0))),
diff --git a/examples/hybrid/types.jl b/examples/hybrid/types.jl
index 9fec96570f..4062734e2d 100644
--- a/examples/hybrid/types.jl
+++ b/examples/hybrid/types.jl
@@ -130,6 +130,8 @@ function get_model_spec(::Type{FT}, parsed_args, namelist) where {FT}
     # should this live in the radiation model?
     idealized_h2o = parsed_args["idealized_h2o"]
     @assert idealized_h2o in (true, false)
+    non_orographic_gravity_wave = parsed_args["non_orographic_gravity_wave"]
+    @assert non_orographic_gravity_wave in (true, false)
 
     model_spec = (;
         moisture_model = moisture_model(parsed_args),
@@ -143,6 +145,7 @@ function get_model_spec(::Type{FT}, parsed_args, namelist) where {FT}
         anelastic_dycore = parsed_args["anelastic_dycore"],
         surface_scheme = surface_scheme(FT, parsed_args),
         C_E = FT(parsed_args["C_E"]),
+        non_orographic_gravity_wave,
     )
 
     return model_spec
diff --git a/perf/flame.jl b/perf/flame.jl
index 27faaa5b84..551f28960a 100644
--- a/perf/flame.jl
+++ b/perf/flame.jl
@@ -67,7 +67,7 @@ allocs = @allocated OrdinaryDiffEq.step!(integrator)
 @info "`allocs ($job_id)`: $(allocs)"
 
 allocs_limit = Dict()
-allocs_limit["flame_perf_target_rhoe"] = 538352
+allocs_limit["flame_perf_target_rhoe"] = 538544
 allocs_limit["flame_perf_target_rhoe_threaded"] = 4440864
 allocs_limit["flame_perf_target_rhoe_callbacks"] = 12832664
 
diff --git a/regression_tests/mse_tables.jl b/regression_tests/mse_tables.jl
index 7041e26254..2337930484 100644
--- a/regression_tests/mse_tables.jl
+++ b/regression_tests/mse_tables.jl
@@ -143,6 +143,13 @@ all_best_mse["compressible_edmf_gabls"][(:c, :uₕ, :components, :data, 1)] = 0.
 all_best_mse["compressible_edmf_gabls"][(:c, :uₕ, :components, :data, 2)] = 0.0
 all_best_mse["compressible_edmf_gabls"][(:c, :turbconv, :en, :ρatke)] = 0.0
 #
+all_best_mse["single_column_nonorographic_gravity_wave"] = OrderedCollections.OrderedDict()
+all_best_mse["single_column_nonorographic_gravity_wave"][(:c, :ρ)] = 0
+all_best_mse["single_column_nonorographic_gravity_wave"][(:c, :ρe_tot)] = 0
+all_best_mse["single_column_nonorographic_gravity_wave"][(:c, :uₕ, :components, :data, 1)] = 0
+all_best_mse["single_column_nonorographic_gravity_wave"][(:c, :uₕ, :components, :data, 2)] = 0
+all_best_mse["single_column_nonorographic_gravity_wave"][(:f, :w, :components, :data, 1)] = 0
+#
 #! format: on
 #################################
 #################################
diff --git a/regression_tests/ref_counter.jl b/regression_tests/ref_counter.jl
index 8351c19397..60d3b2f4a4 100644
--- a/regression_tests/ref_counter.jl
+++ b/regression_tests/ref_counter.jl
@@ -1 +1 @@
-14
+15