From 3df4b53a7c5840962e4e52dbbf9ada6740df035c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E8=B0=A2=E8=90=A7=E6=B6=AF?= Date: Tue, 6 Aug 2024 15:09:04 -0700 Subject: [PATCH] refactor gravity wave and modify tests --- src/cache/temporary_quantities.jl | 2 + .../non_orographic_gravity_wave.jl | 774 ++++++++++++------ .../nogw_test_3d.jl | 101 ++- .../nogw_test_mima.jl | 142 ++-- .../nogw_test_single_column.jl | 231 ++++-- 5 files changed, 876 insertions(+), 374 deletions(-) diff --git a/src/cache/temporary_quantities.jl b/src/cache/temporary_quantities.jl index 04356c41c2..4482e86147 100644 --- a/src/cache/temporary_quantities.jl +++ b/src/cache/temporary_quantities.jl @@ -15,6 +15,8 @@ function temporary_quantities(Y, atmos) ᶜtemp_scalar = Fields.Field(FT, center_space), # ᶜ1 ᶜtemp_scalar_2 = Fields.Field(FT, center_space), # ᶜtke_exch ᶜtemp_scalar_3 = Fields.Field(FT, center_space), + ᶜtemp_scalar_4 = Fields.Field(FT, center_space), + ᶜtemp_scalar_5 = Fields.Field(FT, center_space), ᶠtemp_field_level = Fields.level(Fields.Field(FT, face_space), half), temp_field_level = Fields.level(Fields.Field(FT, center_space), 1), temp_field_level_2 = Fields.level(Fields.Field(FT, center_space), 1), diff --git a/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl b/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl index c0c709baa8..14ee475fd8 100644 --- a/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl +++ b/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl @@ -28,12 +28,13 @@ function non_orographic_gravity_wave_cache( nc = Int(floor(FT(2 * cmax / dc + 1))) c = [FT((n - 1) * dc - cmax) for n in 1:nc] - source_level_z = similar(Fields.level(Y.c.ρ, 1), Tuple{FT, FT}) + source_ρ_z_u_v_level = + similar(Fields.level(Y.c.ρ, 1), Tuple{FT, FT, FT, FT, FT}) ᶜlevel = similar(Y.c.ρ, FT) for i in 1:Spaces.nlevels(axes(Y.c.ρ)) fill!(Fields.level(ᶜlevel, i), i) end - damp_level_z = similar(source_level_z) + return (; gw_source_height = source_height, gw_source_ampl = Bt_0 .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), @@ -41,6 +42,8 @@ function non_orographic_gravity_wave_cache( gw_Bn = Bn .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), gw_B0 = similar(c), gw_c = c, + gw_dc = dc, + gw_cmax = cmax, gw_cw = cw .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), gw_cn = cn .* ones(FT, axes(Fields.level(Y.c.ρ, 1))), gw_c0 = c0, @@ -48,15 +51,15 @@ function non_orographic_gravity_wave_cache( gw_nk = Int(nk), ᶜbuoyancy_frequency = similar(Y.c.ρ), ᶜdTdz = similar(Y.c.ρ), - source_level_z, - damp_level_z, + source_ρ_z_u_v_level, source_level = similar(Fields.level(Y.c.ρ, 1)), damp_level = similar(Fields.level(Y.c.ρ, 1)), ᶜlevel, - u_phy = similar(Y.c.ρ), - v_phy = similar(Y.c.ρ), + u_waveforcing = similar(Y.c.ρ), + v_waveforcing = similar(Y.c.ρ), uforcing = similar(Y.c.ρ), vforcing = similar(Y.c.ρ), + gw_ncval = Val(nc), ) end @@ -82,12 +85,13 @@ function non_orographic_gravity_wave_cache( gw_Bw = ones(FT, axes(lat)) .* Bw gw_cn = ones(FT, axes(lat)) .* cn - source_level_z = similar(Fields.level(Y.c.ρ, 1), Tuple{FT, FT}) + source_p_ρ_z_u_v_level = + similar(Fields.level(Y.c.ρ, 1), Tuple{FT, FT, FT, FT, FT, FT}) ᶜlevel = similar(Y.c.ρ, FT) for i in 1:Spaces.nlevels(axes(Y.c.ρ)) fill!(Fields.level(ᶜlevel, i), i) end - damp_level_z = similar(source_level_z) + # This is GFDL source specs -> a smooth function # source_ampl = @. Bt_0 + @@ -121,20 +125,22 @@ function non_orographic_gravity_wave_cache( gw_c = c, gw_cw = gw_cw, gw_cn = gw_cn, + gw_dc = dc, + gw_cmax = cmax, gw_c0 = c0, gw_flag = gw_flag, gw_nk = Int(nk), ᶜbuoyancy_frequency = similar(Y.c.ρ), ᶜdTdz = similar(Y.c.ρ), - source_level_z, - damp_level_z, + source_p_ρ_z_u_v_level, source_level = similar(Fields.level(Y.c.ρ, 1)), damp_level = similar(Fields.level(Y.c.ρ, 1)), ᶜlevel, - u_phy = similar(Y.c.ρ), - v_phy = similar(Y.c.ρ), + u_waveforcing = similar(Y.c.ρ), + v_waveforcing = similar(Y.c.ρ), uforcing = similar(Y.c.ρ), vforcing = similar(Y.c.ρ), + gw_ncval = Val(nc), ) end @@ -153,33 +159,22 @@ function non_orographic_gravity_wave_tendency!( ᶜdTdz, ᶜbuoyancy_frequency, source_level, - source_level_z, damp_level, - damp_level_z, - u_phy, - v_phy, + u_waveforcing, + v_waveforcing, uforcing, vforcing, ᶜlevel, + gw_ncval, ) = p.non_orographic_gravity_wave (; model_config) = p.atmos - (; - gw_source_ampl, - gw_Bw, - gw_Bn, - gw_B0, - gw_c, - gw_cw, - gw_cn, - gw_flag, - gw_c0, - gw_nk, - ) = p.non_orographic_gravity_wave if model_config isa SingleColumnModel - (; gw_source_height) = p.non_orographic_gravity_wave + (; gw_source_height, source_ρ_z_u_v_level) = + p.non_orographic_gravity_wave elseif model_config isa SphericalModel - (; gw_source_pressure, gw_damp_pressure) = p.non_orographic_gravity_wave + (; gw_source_pressure, gw_damp_pressure, source_p_ρ_z_u_v_level) = + p.non_orographic_gravity_wave end ᶜρ = Y.c.ρ ᶜz = Fields.coordinate_field(Y.c).z @@ -201,23 +196,29 @@ function non_orographic_gravity_wave_tendency!( sqrt(abs(ᶜbuoyancy_frequency)), ) # to avoid small numbers + # prepare physical uv input variables for gravity_wave_forcing() + ᶜu = Geometry.UVVector.(Y.c.uₕ).components.data.:1 + ᶜv = Geometry.UVVector.(Y.c.uₕ).components.data.:2 + if model_config isa SingleColumnModel # source level: the index of the level that is closest to the source height + input = Base.Broadcast.broadcasted(tuple, ᶜρ, ᶜz, ᶜu, ᶜv, ᶜlevel) Operators.column_reduce!( - min_distance_reduce, - source_level_z, - Base.broadcasted( - tuple, - Base.broadcasted( - abs, - Base.broadcasted(-, ᶜz, gw_source_height), - ), - ᶜlevel, - ), - ) + source_ρ_z_u_v_level, + input; + ) do (ρ_prev, z_prev, u_prev, v_prev, level_prev), (ρ, z, u, v, level) + if abs(z_prev - gw_source_height) >= abs(z - gw_source_height) + return (ρ, z, u, v, level) + else + return (ρ_prev, z_prev, u_prev, v_prev, level_prev) + end + end - source_level = source_level_z.:2 + ᶜρ_source = source_ρ_z_u_v_level.:1 + ᶜu_source = source_ρ_z_u_v_level.:3 + ᶜv_source = source_ρ_z_u_v_level.:4 + source_level = source_ρ_z_u_v_level.:5 fill!(damp_level, Spaces.nlevels(axes(ᶜz))) @@ -225,228 +226,511 @@ function non_orographic_gravity_wave_tendency!( (; ᶜp) = p.precomputed # source level: the index of the highest level whose pressure is higher than source pressure + input = Base.Broadcast.broadcasted(tuple, ᶜp, ᶜρ, ᶜz, ᶜu, ᶜv, ᶜlevel) Operators.column_reduce!( - positive_selector_reduce, - source_level_z, - Base.broadcasted( - tuple, - Base.broadcasted(-, ᶜp, gw_source_pressure), - ᶜlevel, - ), - ) + source_p_ρ_z_u_v_level, + input, + ) do (p_prev, ρ_prev, z_prev, u_prev, v_prev, level_prev), + (p, ρ, z, u, v, level) + if (p - gw_source_pressure) <= 0 + return (p_prev, ρ_prev, z_prev, u_prev, v_prev, level_prev) + else + return (p, ρ, z, u, v, level) + end + end - source_level = source_level_z.:2 + ᶜρ_source = source_p_ρ_z_u_v_level.:2 + ᶜu_source = source_p_ρ_z_u_v_level.:4 + ᶜv_source = source_p_ρ_z_u_v_level.:5 + source_level = source_p_ρ_z_u_v_level.:6 # damp level: the index of the lowest level whose pressure is lower than the damp pressure + input = Base.Broadcast.broadcasted(tuple, ᶜlevel, ᶜp) Operators.column_reduce!( - negative_selector_reduce, - damp_level_z, - Base.broadcasted( - tuple, - Base.broadcasted(-, ᶜp, gw_damp_pressure), - ᶜlevel, - ), - ) - - damp_level = damp_level_z.:2 + damp_level, + input; + transform = first, + ) do (level_prev, p_prev), (level, p) + if (p_prev - gw_damp_pressure) >= 0 + return (level, p) + else + return (level_prev, p_prev) + end + end end - # prepare physical uv input variables for gravity_wave_forcing() - u_phy = Geometry.UVVector.(Y.c.uₕ).components.data.:1 - v_phy = Geometry.UVVector.(Y.c.uₕ).components.data.:2 - - # GW parameterization applied bycolume - Fields.bycolumn(axes(ᶜρ)) do colidx - parent(uforcing[colidx]) .= non_orographic_gravity_wave_forcing( - vec(parent(u_phy[colidx])), - vec(parent(ᶜbuoyancy_frequency[colidx])), - vec(parent(ᶜρ[colidx])), - vec(parent(ᶜz[colidx])), - Int(parent(source_level[colidx])[1]), - Int(parent(damp_level[colidx])[1]), - parent(gw_source_ampl[colidx])[1], - parent(gw_Bw[colidx])[1], - parent(gw_Bn[colidx])[1], - gw_B0, - parent(gw_cw[colidx])[1], - parent(gw_cn[colidx])[1], - parent(gw_flag[colidx])[1], - gw_c, - gw_c0, - gw_nk, - ) + ᶜu = Geometry.UVVector.(Y.c.uₕ).components.data.:1 + ᶜv = Geometry.UVVector.(Y.c.uₕ).components.data.:2 - parent(vforcing[colidx]) .= non_orographic_gravity_wave_forcing( - vec(parent(v_phy[colidx])), - vec(parent(ᶜbuoyancy_frequency[colidx])), - vec(parent(ᶜρ[colidx])), - vec(parent(ᶜz[colidx])), - Int(parent(source_level[colidx])[1]), - Int(parent(damp_level[colidx])[1]), - parent(gw_source_ampl[colidx])[1], - parent(gw_Bw)[1], - parent(gw_Bn)[1], - gw_B0, - parent(gw_cw)[1], - parent(gw_cn)[1], - parent(gw_flag)[1], - gw_c, - gw_c0, - gw_nk, - ) + uforcing .= 0 + vforcing .= 0 - end + non_orographic_gravity_wave_forcing( + ᶜu, + ᶜv, + ᶜbuoyancy_frequency, + ᶜρ, + ᶜz, + ᶜlevel, + source_level, + damp_level, + ᶜρ_source, + ᶜu_source, + ᶜv_source, + uforcing, + vforcing, + gw_ncval, + u_waveforcing, + v_waveforcing, + p, + ) - # physical uv forcing converted to Covariant12Vector and added up to uₕ tendencies @. Yₜ.c.uₕ += Geometry.Covariant12Vector.(Geometry.UVVector.(uforcing, vforcing)) - return nothing + end function non_orographic_gravity_wave_forcing( - old_ᶜu, - old_ᶜbf, - old_ᶜρ, - old_ᶜz, + ᶜu, + ᶜv, + ᶜbf, + ᶜρ, + ᶜz, + ᶜlevel, source_level, damp_level, - source_ampl, - Bw, - Bn, - B0, - cw, - cn, - flag, - c, - c0, - nk, -) - FT = eltype(old_ᶜz) - # add an extra layer above model top so that forcing between the very top - # model layer and the upper boundary can be calculated - ᶜu = vcat(old_ᶜu, FT(2) * old_ᶜu[end] - old_ᶜu[end - 1]) - ᶜρ = vcat(old_ᶜρ, old_ᶜρ[end] * old_ᶜρ[end] / old_ᶜρ[end - 1]) - ᶜbf = vcat(old_ᶜbf, old_ᶜbf[end]) - ᶜz = vcat(old_ᶜz, FT(2) * old_ᶜz[end] - old_ᶜz[end - 1]) - - # wave spectra and the source amplitude - nc = length(c) - c_hat0 = c .- ᶜu[source_level] # c0mu0 - - Bw_exp = @. exp(-log(2.0) * ((c * flag + c_hat0 * (1 - flag) - c0) / cw)^2) - Bn_exp = @. exp(-log(2.0) * ((c * flag + c_hat0 * (1 - flag) - c0) / cn)^2) - B0 = @. sign(c_hat0) * (Bw * Bw_exp + Bn * Bn_exp) - - Bsum = sum(abs.(B0)) - if (Bsum == 0.0) - error("zero flux input at source level") - end - # intermittency - eps = calc_intermitency(ᶜρ[source_level], source_ampl, nk, Bsum) - - # horizontal wave length - kwv = [2.0 * π / ((30.0 * (10.0^n)) * 1.e3) for n in 1:nk] - k2 = kwv .* kwv - - # forcing - wave_forcing = zeros(length(ᶜu)) - gwf = zeros(length(ᶜu) - 1) - 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) # here ᶜu has one additional level above model top - fac = FT(0.5) * (ᶜρ[k] / ᶜρ[source_level]) * kwv[ink] / ᶜbf[k] - - ᶜHb = -(ᶜz[k] - ᶜz[k - 1]) / log(ᶜρ[k] / ᶜρ[k - 1]) # density scale height - alp2 = 0.25 / (ᶜHb * ᶜHb) - ω_r = sqrt((ᶜbf[k] * ᶜbf[k] * k2[ink]) / (k2[ink] + alp2)) # omc: (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] # c0mu - # 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 + ᶜρ_source, + ᶜu_source, + ᶜv_source, + uforcing, + vforcing, + gw_ncval::Val{nc}, + u_waveforcing, + v_waveforcing, + p, +) where {nc} + (; + gw_source_ampl, + gw_Bw, + gw_Bn, + gw_c, + gw_cw, + gw_cn, + gw_flag, + gw_c0, + gw_nk, + ) = p.non_orographic_gravity_wave + ᶜρ_p1 = p.scratch.ᶜtemp_scalar + ᶜz_p1 = p.scratch.ᶜtemp_scalar_2 + ᶜu_p1 = p.scratch.ᶜtemp_scalar_3 + ᶜv_p1 = p.scratch.ᶜtemp_scalar_4 + ᶜbf_p1 = p.scratch.ᶜtemp_scalar_5 + nci = get_nc(gw_ncval) + FT = eltype(ᶜρ) + ρ_endlevel = Fields.level(ᶜρ, Spaces.nlevels(axes(ᶜρ))) + ρ_endlevel_m1 = Fields.level(ᶜρ, Spaces.nlevels(axes(ᶜρ)) - 1) + Boundary_value = Fields.Field( + Fields.field_values(ρ_endlevel) .* Fields.field_values(ρ_endlevel) ./ + Fields.field_values(ρ_endlevel_m1), + axes(ρ_endlevel), + ) + field_shiftlevel_up!(ᶜρ, ᶜρ_p1, Boundary_value) + + u_endlevel = Fields.level(ᶜu, Spaces.nlevels(axes(ᶜu))) + u_endlevel_m1 = Fields.level(ᶜu, Spaces.nlevels(axes(ᶜu)) - 1) + Boundary_value = Fields.Field( + FT(2) .* Fields.field_values(u_endlevel) .- + Fields.field_values(u_endlevel_m1), + axes(u_endlevel), + ) + field_shiftlevel_up!(ᶜu, ᶜu_p1, Boundary_value) + + v_endlevel = Fields.level(ᶜv, Spaces.nlevels(axes(ᶜv))) + v_endlevel_m1 = Fields.level(ᶜv, Spaces.nlevels(axes(ᶜv)) - 1) + Boundary_value = Fields.Field( + FT(2) .* Fields.field_values(v_endlevel) .- + Fields.field_values(v_endlevel_m1), + axes(v_endlevel), + ) + field_shiftlevel_up!(ᶜv, ᶜv_p1, Boundary_value) + + Boundary_value = Fields.level(ᶜbf, Spaces.nlevels(axes(ᶜbf))) + field_shiftlevel_up!(ᶜbf, ᶜbf_p1, Boundary_value) + + z_endlevel = Fields.level(ᶜz, Spaces.nlevels(axes(ᶜz))) + z_endlevel_m1 = Fields.level(ᶜz, Spaces.nlevels(axes(ᶜz)) - 1) + Boundary_value = Fields.Field( + FT(2) .* Fields.field_values(z_endlevel) .- + Fields.field_values(z_endlevel_m1), + axes(z_endlevel), + ) + field_shiftlevel_up!(ᶜz, ᶜz_p1, Boundary_value) + + mask = BitVector(ones(nc)) + level_end = Spaces.nlevels(axes(ᶜρ)) + B1 = ntuple(i -> 0.0, Val(nc)) + + for ink in 1:gw_nk + + input_u = Base.Broadcast.broadcasted( + tuple, + ᶜu_p1, + ᶜu_source, + ᶜbf_p1, + ᶜρ, + ᶜρ_p1, + ᶜρ_source, + ᶜz_p1, + ᶜz, + source_level, + gw_Bw, + gw_Bn, + gw_cw, + gw_cn, + gw_c, + gw_flag, + ᶜlevel, + gw_source_ampl, + ) + + + input_v = Base.Broadcast.broadcasted( + tuple, + ᶜv_p1, + ᶜv_source, + ᶜbf_p1, + ᶜρ, + ᶜρ_p1, + ᶜρ_source, + ᶜz_p1, + ᶜz, + source_level, + gw_Bw, + gw_Bn, + gw_cw, + gw_cn, + gw_c, + gw_flag, + ᶜlevel, + gw_source_ampl, + ) + + Operators.column_accumulate!( + u_waveforcing, + input_u; + init = (FT(0.0), mask, 0.0, B1), + transform = first, + ) do (wave_forcing, mask, Bsum, B0), + ( + u_kp1, + u_source, + bf_kp1, + ρ_k, + ρ_kp1, + ρ_source, + z_kp1, + z_k, + source_level, + Bw, + Bn, + cw, + cn, + c, + flag, + level, + source_ampl, + ) + if level >= (source_level - 1) + FT1 = typeof(u_kp1) + kwv = 2.0 * π / ((30.0 * (10.0^ink)) * 1.e3) + k2 = kwv * kwv + fac = FT1(0.5) * (ρ_kp1 / ρ_source) * kwv / bf_kp1 + Hb = (z_kp1 - z_k) / log(ρ_k / ρ_kp1) # density scale height + alp2 = FT1(0.25) / (Hb * Hb) + ω_r = sqrt((bf_kp1 * bf_kp1 * k2) / (k2 + alp2)) # omc: (critical frequency that marks total internal reflection) + fm = 0.0 + if level == (source_level - 1) + mask .= 1 + Bsum = 0.0 + B0 = ntuple( + n -> + sign(c[n] - u_source) * ( + Bw * exp( + -log(2.0) * + ( + ( + c[n] * flag + + (c[n] - u_source) * (1 - flag) - + gw_c0 + ) / cw + )^2, + ) + + Bn * exp( + -log(2.0) * + ( + ( + c[n] * flag + + (c[n] - u_source) * (1 - flag) - + gw_c0 + ) / cn + )^2, + ) + ), + Val(nc), + ) + Bsum = sum(abs.(B0)) + end + for n in 1:nci + # check only those waves which are still propagating, i.e., mask = 1.0 + if (mask[n]) == 1 + c_hat = c[n] - u_kp1 # c0mu + # f phase speed matches the wind speed, remove c(n) from the set of propagating waves. + if c_hat == 0.0 + mask[n] = 0 else - if k == length(ᶜu) - # this is added in MiMA implementation: - # all momentum flux that escapes across the model top - # is deposited to the extra level being added so that - # momentum flux is conserved - mask[n] = 0.0 - if k > source_level - fm = fm + B0[n] - end + c_hat0 = c[n] - u_source + # define the criterion which determines if wave is reflected at this level (test). + test = abs(c_hat) * kwv - ω_r + if test >= 0.0 + # wave has undergone total internal reflection. remove it from the propagating set. + mask[n] = 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. - Foc = B0[n] / (c_hat)^3 - fac - if Foc >= 0.0 || (c_hat0[n] * c_hat <= 0.0) - mask[n] = 0.0 - if k > source_level + if level == level_end + # this is added in MiMA implementation: + # all momentum flux that escapes across the model top + # is deposited to the extra level being added so that + # momentum flux is conserved + mask[n] = 0 + if level >= source_level fm = fm + B0[n] end + 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. + Foc = B0[n] / FT1((c_hat)^3) - fac + if Foc >= 0.0 || (c_hat0 * c_hat <= 0.0) + mask[n] = 0 + if level >= source_level + fm = fm + B0[n] + end + end end - end - end # (test >= 0.0) - end #(c_hat == 0.0) - end # mask = 0 - - end # nc: phase speed loop - - # 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 / (ᶜz[k] - ᶜz[k - 1]) - if k == length(ᶜu) - wave_forcing[k - 1] = 0.5 * wave_forcing[k - 1] + end # (test >= 0.0) + end #(c_hat == 0.0) + end # mask = 0 + + end # nc: phase speed loop + + # compute the gravity wave momentum flux forcing + # obtained across the entire wave spectrum at this level. + eps = calc_intermitency(ρ_source, source_ampl, gw_nk, FT1(Bsum)) + if level >= source_level + rbh = sqrt(ρ_k * ρ_kp1) + wave_forcing = + (ρ_source / rbh) * FT1(fm) * eps / (z_kp1 - z_k) else - wave_forcing[k - 1] = - 0.5 * (wave_forcing[k - 1] + wave_forcing[k]) + wave_forcing = FT1(0.0) end - else - wave_forcing[k] = 0.0 end + return (wave_forcing, mask, Bsum, B0) - end # k - - # model top: deposit remaining momentum flux that goes across the model top - # to the levels above the damp level - # This is not included in Joan Alexander's original code nor the GFDL implementation; - # but is added in MiMA based on Tiffany Shaw's paper: - # https://journals.ametsoc.org/view/journals/clim/22/10/2009jcli2688.1.xml?tab_body=pdf - for k in damp_level:(length(ᶜu) - 1) - wave_forcing[k] = - wave_forcing[k] + - wave_forcing[end] / (length(ᶜu) + 1 - damp_level) end - # forcing - for k in source_level:(length(ᶜu) - 1) - gwf[k] = gwf[k] + wave_forcing[k] + + + Operators.column_accumulate!( + v_waveforcing, + input_v; + init = (FT(0.0), mask, 0.0, B1), + transform = first, + ) do (wave_forcing, mask, Bsum, B0), + ( + u_kp1, + u_source, + bf_kp1, + ρ_k, + ρ_kp1, + ρ_source, + z_kp1, + z_k, + source_level, + Bw, + Bn, + cw, + cn, + c, + flag, + level, + source_ampl, + ) + + if level >= (source_level - 1) + FT2 = typeof(u_kp1) + kwv = 2.0 * π / ((30.0 * (10.0^ink)) * 1.e3) + k2 = kwv * kwv + fac = FT2(0.5) * (ρ_kp1 / ρ_source) * kwv / bf_kp1 + Hb = (z_kp1 - z_k) / log(ρ_k / ρ_kp1) # density scale height + alp2 = FT2(0.25) / (Hb * Hb) + ω_r = sqrt((bf_kp1 * bf_kp1 * k2) / (k2 + alp2)) # omc: (critical frequency that marks total internal reflection) + + fm = 0.0 + if level == (source_level - 1) + mask .= 1 + Bsum = 0.0 + B0 = ntuple( + n -> + sign((c[n] - u_source)) * ( + Bw * exp( + -log(2.0) * + ( + ( + c[n] * flag + + (c[n] - u_source) * (1 - flag) - + gw_c0 + ) / cw + )^2, + ) + + Bn * exp( + -log(2.0) * + ( + ( + c[n] * flag + + (c[n] - u_source) * (1 - flag) - + gw_c0 + ) / cn + )^2, + ) + ), + Val(nc), + ) + Bsum = sum(abs.(B0)) + end + for n in 1:nci + # check only those waves which are still propagating, i.e., mask = 1.0 + if (mask[n]) == 1 + c_hat = c[n] - u_kp1 # c0mu + # f phase speed matches the wind speed, remove c(n) from the set of propagating waves. + if c_hat == 0.0 + mask[n] = 0 + else + c_hat0 = c[n] - u_source + # define the criterion which determines if wave is reflected at this level (test). + test = abs(c_hat) * kwv - ω_r + if test >= 0.0 + # wave has undergone total internal reflection. remove it from the propagating set. + mask[n] = 0 + else + if level == level_end + # this is added in MiMA implementation: + # all momentum flux that escapes across the model top + # is deposited to the extra level being added so that + # momentum flux is conserved + mask[n] = 0 + if level >= source_level + fm = fm + B0[n] + end + 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. + Foc = B0[n] / FT2((c_hat)^3) - fac + if Foc >= 0.0 || (c_hat0 * c_hat <= 0.0) + mask[n] = 0 + if level >= source_level + fm = fm + B0[n] + end + end + end + end # (test >= 0.0) + end #(c_hat == 0.0) + end # mask = 0 + + end # nc: phase speed loop + + # compute the gravity wave momentum flux forcing + # obtained across the entire wave spectrum at this level. + eps = calc_intermitency(ρ_source, source_ampl, gw_nk, FT2(Bsum)) + if level >= source_level + rbh = sqrt(ρ_k * ρ_kp1) + wave_forcing = + (ρ_source / rbh) * FT2(fm) * eps / (z_kp1 - z_k) + else + wave_forcing = FT2(0.0) + end + end + return (wave_forcing, mask, Bsum, B0) + end - end # ink + u_waveforcing_top = p.scratch.temp_field_level + copyto!( + Fields.field_values(u_waveforcing_top), + Fields.field_values( + Fields.level( + u_waveforcing, + Spaces.nlevels(axes(u_waveforcing)), + ), + ), + ) + fill!( + Fields.level(u_waveforcing, Spaces.nlevels(axes(u_waveforcing))), + 0, + ) - return gwf + # v_waveforcing_top = similar( + # Fields.level(v_waveforcing, Spaces.nlevels(axes(v_waveforcing))), + # ) + v_waveforcing_top = p.scratch.temp_field_level + copyto!( + Fields.field_values(v_waveforcing_top), + Fields.field_values( + Fields.level( + v_waveforcing, + Spaces.nlevels(axes(v_waveforcing)), + ), + ), + ) + fill!( + Fields.level(v_waveforcing, Spaces.nlevels(axes(v_waveforcing))), + 0, + ) + + gw_average!(u_waveforcing, p.scratch.ᶜtemp_scalar) + gw_average!(v_waveforcing, p.scratch.ᶜtemp_scalar) + + @. u_waveforcing = gw_deposit( + u_waveforcing_top, + u_waveforcing, + damp_level, + ᶜlevel, + level_end, + ) + @. v_waveforcing = gw_deposit( + v_waveforcing_top, + v_waveforcing, + damp_level, + ᶜlevel, + level_end, + ) + + @. uforcing = uforcing + u_waveforcing + @. vforcing = vforcing + v_waveforcing + + end + return nothing end # calculate the intermittency factor eps -> assuming constant Δc. @@ -455,6 +739,28 @@ function calc_intermitency(ρ_source_level, source_ampl, nk, Bsum) return (source_ampl / ρ_source_level / nk) / Bsum end -@inline min_distance_reduce(a, b) = ifelse(a[1] < b[1], a, b) -@inline positive_selector_reduce(a, b) = ifelse(b[1] <= 0, a, b) -@inline negative_selector_reduce(a, b) = ifelse(a[1] >= 0, b, a) +function gw_average!(wave_forcing, wave_forcing_m1) + L1 = Operators.LeftBiasedC2F(; bottom = Operators.SetValue(0.0)) + L2 = Operators.LeftBiasedF2C(;) + wave_forcing_m1 .= L2.(L1.(wave_forcing)) + @. wave_forcing = 0.5 * (wave_forcing + wave_forcing_m1) +end + + +function gw_deposit(wave_forcing_top, wave_forcing, damp_level, level, height) + if level >= damp_level + wave_forcing = + wave_forcing + wave_forcing_top / (height + 2 - damp_level) + end + return wave_forcing +end + +function get_nc(::Val{nc}) where {nc} + nc +end + +function field_shiftlevel_up!(ᶜexample_field, ᶜshifted_field, Boundary_value) + R1 = Operators.RightBiasedC2F(; top = Operators.SetValue(Boundary_value)) + R2 = Operators.RightBiasedF2C(;) + ᶜshifted_field .= R2.(R1.(ᶜexample_field)) +end diff --git a/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl b/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl index 45cbe75b0c..42d3829e12 100644 --- a/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl +++ b/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl @@ -6,6 +6,11 @@ using Interpolations using Statistics import ClimaAtmos import ClimaAtmos as CA +import ClimaCore +import ClimaCore.Spaces as Spaces +import ClimaCore.Fields as Fields +import ClimaCore.Geometry as Geometry +import ClimaCore.Operators as Operators const FT = Float64 include("../gw_plotutils.jl") @@ -17,7 +22,7 @@ face_z = FT.(0:1e3:0.47e5) center_z = FT(0.5) .* (face_z[1:(end - 1)] .+ face_z[2:end]) # compute the source parameters -function non_orographic_gravity_wave( +function non_orographic_gravity_wave_param( ::Type{FT}; source_height = FT(15000), Bw = FT(1.2), @@ -46,8 +51,14 @@ function non_orographic_gravity_wave( ) end -params = non_orographic_gravity_wave(FT; Bw = 0.4, cmax = 150, kwv = 2π / 100e3) -source_level = argmin(abs.(center_z .- params.gw_source_height)) +non_orographic_gravity_wave = non_orographic_gravity_wave_param( + FT; + Bw = 0.4, + cmax = 150, + kwv = 2π / 100e3, +) +source_level = + argmin(abs.(center_z .- non_orographic_gravity_wave.gw_source_height)) damp_level = length(center_z) include(joinpath(pkgdir(ClimaAtmos), "artifacts", "artifact_funcs.jl")) @@ -122,7 +133,45 @@ 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, :, :, :] -B0 = similar(params.gw_c) +column_domain = ClimaCore.Domains.IntervalDomain( + ClimaCore.Geometry.ZPoint(0.0) .. ClimaCore.Geometry.ZPoint(47000), + boundary_names = (:bottom, :top), +) + +column_mesh = ClimaCore.Meshes.IntervalMesh(column_domain, nelems = 47) + +column_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(column_mesh) + +column_face_space = + ClimaCore.Spaces.FaceFiniteDifferenceSpace(column_center_space) + +coordinate = Fields.coordinate_field(column_center_space) +gw_ncval = Val(500) +ᶜz = coordinate.z +ᶜρ = copy(ᶜz) +ᶜu = copy(ᶜz) +ᶜv = copy(ᶜz) +ᶜbf = copy(ᶜz) +ᶜlevel = similar(ᶜρ, FT) +u_waveforcing = similar(ᶜu) +v_waveforcing = similar(ᶜv) +for i in 1:Spaces.nlevels(axes(ᶜρ)) + fill!(Fields.level(ᶜlevel, i), i) +end +uforcing = similar(ᶜρ, FT) +vforcing = similar(ᶜρ, FT) + +scratch = (; + ᶜtemp_scalar = similar(ᶜz, FT), + ᶜtemp_scalar_2 = similar(ᶜz, FT), + ᶜtemp_scalar_3 = similar(ᶜz, FT), + ᶜtemp_scalar_4 = similar(ᶜz, FT), + ᶜtemp_scalar_5 = similar(ᶜz, FT), + temp_field_level = similar(Fields.level(ᶜz, 1), FT), +) + +params = (; non_orographic_gravity_wave, scratch) + # Jan month = Dates.month.(time) @@ -130,26 +179,38 @@ 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)) +Jan_uforcing = similar(Jan_u) + for j in 1:length(lat) - Jan_uforcing[j, :] = CA.non_orographic_gravity_wave_forcing( - Jan_u[j, :], - Jan_bf[j, :], - Jan_ρ[j, :], - copy(center_z), + Base.parent(ᶜρ) .= Jan_ρ[j, :] + Base.parent(ᶜu) .= Jan_u[j, :] + Base.parent(ᶜbf) .= Jan_bf[j, :] + ᶜρ_source = Fields.level(ᶜρ, source_level) + ᶜu_source = Fields.level(ᶜu, source_level) + ᶜv_source = Fields.level(ᶜv, source_level) + uforcing .= 0 + vforcing .= 0 + + CA.non_orographic_gravity_wave_forcing( + ᶜu, + ᶜv, + ᶜbf, + ᶜρ, + ᶜz, + ᶜlevel, source_level, damp_level, - params.gw_source_ampl, - params.gw_Bw, - params.gw_Bn, - B0, - params.gw_cw, - params.gw_cn, - params.gw_flag, - params.gw_c, - params.gw_c0, - params.gw_nk, + ᶜρ_source, + ᶜu_source, + ᶜv_source, + uforcing, + vforcing, + gw_ncval, + u_waveforcing, + v_waveforcing, + params, ) + Jan_uforcing[j, :] = parent(uforcing) end output_dir = "nonorographic_gravity_wave_test_3d" diff --git a/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl b/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl index 13719a2589..d44f92c09b 100644 --- a/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl +++ b/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl @@ -6,12 +6,17 @@ using Interpolations using Statistics import ClimaAtmos import ClimaAtmos as CA +import ClimaCore +import ClimaCore.Spaces as Spaces +import ClimaCore.Fields as Fields +import ClimaCore.Geometry as Geometry +import ClimaCore.Operators as Operators const FT = Float64 include("../gw_plotutils.jl") # compute the source parameters -function non_orographic_gravity_wave( +function non_orographic_gravity_wave_param( lat, ::Type{FT}; gw_source_pressure = FT(31500), @@ -75,7 +80,7 @@ function non_orographic_gravity_wave( gw_cw = gw_cw, gw_cn = cn, gw_c0 = c0, - gw_flag = gw_flag, + gw_flag = 0, gw_nk = nk, ) end @@ -124,53 +129,63 @@ bf = @. (grav / T) * (dTdz + grav / cp_d) bf = @. ifelse(bf < 2.5e-5, sqrt(2.5e-5), sqrt(abs(bf))) # compute u/v forcings from convective gravity waves -params = non_orographic_gravity_wave(lat, FT) -B0 = similar(params.gw_c) +param = non_orographic_gravity_wave_param(lat, FT) # nogw forcing kmax = length(pfull) - 1 -k_source = findfirst(pfull * 100 .> params.gw_source_pressure) -k_damp = findlast(pfull * 100 .< params.gw_damp_pressure) +k_source = findfirst(pfull * 100 .> param.gw_source_pressure) +k_damp = findlast(pfull * 100 .< param.gw_damp_pressure) uforcing = zeros(size(lev)) vforcing = zeros(size(lev)) -for i in 1:length(lon) - for j in 1:length(lat) - for it in 1:length(time) - source_level = - kmax + 2 - Int( - floor( - (kmax + 1) - - ((kmax + 1 - k_source) * cosd(lat[j]) + 0.5), - ), - ) - damp_level = length(pfull) + 1 - k_damp - uforcing[i, j, :, it] = CA.non_orographic_gravity_wave_forcing( - u[i, j, end:-1:1, it], - bf[i, j, end:-1:1, it], - ρ[i, j, end:-1:1, it], - z[i, j, end:-1:1, it], - source_level, - damp_level, - params.gw_source_ampl[j], - params.gw_Bw, - params.gw_Bn[j], - B0, - params.gw_cw[j], - params.gw_cn, - params.gw_flag[j], - params.gw_c, - params.gw_c0, - params.gw_nk, - ) - end - end +column_domain = ClimaCore.Domains.IntervalDomain( + ClimaCore.Geometry.ZPoint(FT(z[1, 1, end, 1])) .. + ClimaCore.Geometry.ZPoint(FT(z[1, 1, 1, 1])), + boundary_names = (:bottom, :top), +) + +column_mesh = ClimaCore.Meshes.IntervalMesh(column_domain, nelems = 40) + +# construct the face space from the center one +column_face_space = ClimaCore.Spaces.FaceFiniteDifferenceSpace(column_mesh) +column_center_space = + ClimaCore.Spaces.CenterFiniteDifferenceSpace(column_face_space) + +coordinate = Fields.coordinate_field(column_center_space) + +coordinate = Fields.coordinate_field(column_center_space) +gw_ncval = Val(333) +ᶜz = coordinate.z +ᶜρ = copy(ᶜz) +ᶜu = copy(ᶜz) +ᶜv = copy(ᶜz) +ᶜbf = copy(ᶜz) +ᶜlevel = similar(ᶜρ, FT) +u_waveforcing = similar(ᶜu) +v_waveforcing = similar(ᶜv) +for i in 1:Spaces.nlevels(axes(ᶜρ)) + fill!(Fields.level(ᶜlevel, i), i) end -uforcing_zonalave = dropdims(mean(uforcing, dims = 1), dims = 1) +ᶜuforcing = similar(ᶜρ, FT) +ᶜvforcing = similar(ᶜρ, FT) + -for i in 1:length(lon) - for j in 1:length(lat) +scratch = (; + ᶜtemp_scalar = similar(ᶜz, FT), + ᶜtemp_scalar_2 = similar(ᶜz, FT), + ᶜtemp_scalar_3 = similar(ᶜz, FT), + ᶜtemp_scalar_4 = similar(ᶜz, FT), + ᶜtemp_scalar_5 = similar(ᶜz, FT), + temp_field_level = similar(Fields.level(ᶜz, 1), FT), +) + +#params = (; non_orographic_gravity_wave, scratch) + +for j in 1:length(lat) + non_orographic_gravity_wave = non_orographic_gravity_wave_param(lat[j], FT) + params = (; non_orographic_gravity_wave, scratch) + for i in 1:length(lon) for it in 1:length(time) source_level = kmax + 2 - Int( @@ -180,28 +195,43 @@ for i in 1:length(lon) ), ) damp_level = length(pfull) + 1 - k_damp - - vforcing[i, j, :, it] = CA.non_orographic_gravity_wave_forcing( - v[i, j, end:-1:1, it], - bf[i, j, end:-1:1, it], - ρ[i, j, end:-1:1, it], - z[i, j, end:-1:1, it], + Base.parent(ᶜz) .= z[i, j, end:-1:1, it] + Base.parent(ᶜρ) .= ρ[i, j, end:-1:1, it] + Base.parent(ᶜu) .= u[i, j, end:-1:1, it] + Base.parent(ᶜv) .= v[i, j, end:-1:1, it] + Base.parent(ᶜbf) .= bf[i, j, end:-1:1, it] + ᶜρ_source = Fields.level(ᶜρ, source_level) + ᶜu_source = Fields.level(ᶜu, source_level) + ᶜv_source = Fields.level(ᶜv, source_level) + ᶜuforcing .= 0 + ᶜvforcing .= 0 + + CA.non_orographic_gravity_wave_forcing( + ᶜu, + ᶜv, + ᶜbf, + ᶜρ, + ᶜz, + ᶜlevel, source_level, damp_level, - params.gw_source_ampl[j], - params.gw_Bw, - params.gw_Bn[j], - B0, - params.gw_cw[j], - params.gw_cn, - params.gw_flag[j], - params.gw_c, - params.gw_c0, - params.gw_nk, + ᶜρ_source, + ᶜu_source, + ᶜv_source, + ᶜuforcing, + ᶜvforcing, + gw_ncval, + u_waveforcing, + v_waveforcing, + params, ) + + uforcing[i, j, :, it] = parent(ᶜuforcing) + vforcing[i, j, :, it] = parent(ᶜvforcing) end end end +uforcing_zonalave = dropdims(mean(uforcing, dims = 1), dims = 1) vforcing_zonalave = dropdims(mean(vforcing, dims = 1), dims = 1) # plots diff --git a/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl b/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl index aa98a18a51..e4aab21a26 100644 --- a/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl +++ b/test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl @@ -6,6 +6,11 @@ using Interpolations using Statistics import ClimaAtmos import ClimaAtmos as CA +import ClimaCore +import ClimaCore.Spaces as Spaces +import ClimaCore.Fields as Fields +import ClimaCore.Geometry as Geometry +import ClimaCore.Operators as Operators include("../gw_plotutils.jl") @@ -18,7 +23,7 @@ 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 non_orographic_gravity_wave( +function non_orographic_gravity_wave_param( ::Type{FT}; source_height = FT(15000), Bw = FT(1.2), @@ -47,8 +52,14 @@ function non_orographic_gravity_wave( ) end -params = non_orographic_gravity_wave(FT; Bw = 0.4, cmax = 150, kwv = 2π / 100e3) -source_level = argmin(abs.(center_z .- params.gw_source_height)) +non_orographic_gravity_wave = non_orographic_gravity_wave_param( + FT; + Bw = 0.4, + cmax = 150, + kwv = 2π / 100e3, +) +source_level = + argmin(abs.(center_z .- non_orographic_gravity_wave.gw_source_height)) damp_level = length(center_z) include(joinpath(pkgdir(ClimaAtmos), "artifacts", "artifact_funcs.jl")) @@ -122,6 +133,40 @@ 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 + +column_domain = ClimaCore.Domains.IntervalDomain( + ClimaCore.Geometry.ZPoint(0.0) .. ClimaCore.Geometry.ZPoint(50000.0), + boundary_names = (:bottom, :top), +) + +column_mesh = ClimaCore.Meshes.IntervalMesh(column_domain, nelems = 50) + +column_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(column_mesh) +# construct the face space from the center one +column_face_space = + ClimaCore.Spaces.FaceFiniteDifferenceSpace(column_center_space) + +coord = ClimaCore.Fields.coordinate_field(column_center_space) + +gw_ncval = Val(501) +ᶜz = coord.z +ᶜρ = copy(ᶜz) +ᶜu = copy(ᶜz) +ᶜv = copy(ᶜz) +ᶜbf = copy(ᶜz) +ᶜlevel = similar(ᶜρ, FT) +u_waveforcing = similar(ᶜu) +v_waveforcing = similar(ᶜv) +for i in 1:Spaces.nlevels(axes(ᶜρ)) + fill!(Fields.level(ᶜlevel, i), i) +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) @@ -129,30 +174,53 @@ ENV["GKSwstype"] = "nul" output_dir = "nonorographic_gravity_wave_test_single_column" mkpath(output_dir) -B0 = similar(params.gw_c) +scratch = (; + ᶜtemp_scalar = similar(ᶜz, FT), + ᶜtemp_scalar_2 = similar(ᶜz, FT), + ᶜtemp_scalar_3 = similar(ᶜz, FT), + ᶜtemp_scalar_4 = similar(ᶜz, FT), + ᶜtemp_scalar_5 = similar(ᶜz, FT), + temp_field_level = similar(Fields.level(ᶜz, 1), FT), +) + +params = (; non_orographic_gravity_wave, scratch) # 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 = CA.non_orographic_gravity_wave_forcing( - Jan_u, - Jan_bf, - Jan_ρ, - center_z, +Base.parent(ᶜρ) .= Jan_ρ +ᶜv = ᶜu +Base.parent(ᶜu) .= Jan_u +Base.parent(ᶜbf) .= Jan_bf +ᶜρ_source = Fields.level(ᶜρ, source_level) +ᶜu_source = Fields.level(ᶜu, source_level) +ᶜv_source = Fields.level(ᶜv, source_level) +Jan_uforcing = similar(ᶜρ, FT) +Jan_uforcing .= 0 +Jan_vforcing = similar(ᶜρ, FT) +Jan_vforcing .= 0 + +CA.non_orographic_gravity_wave_forcing( + ᶜu, + ᶜv, + ᶜbf, + ᶜρ, + ᶜz, + ᶜlevel, source_level, damp_level, - params.gw_source_ampl, - params.gw_Bw, - params.gw_Bn, - B0, - params.gw_cw, - params.gw_cn, - params.gw_flag, - params.gw_c, - params.gw_c0, - params.gw_nk, + ᶜρ_source, + ᶜu_source, + ᶜv_source, + Jan_uforcing, + Jan_vforcing, + gw_ncval, + u_waveforcing, + v_waveforcing, + params, ) +Jan_uforcing = Base.parent(Jan_uforcing) fig = generate_empty_figure(); create_plot!( fig; @@ -166,24 +234,35 @@ CairoMakie.save(joinpath(output_dir, "fig6jan.png"), fig) 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 = CA.non_orographic_gravity_wave_forcing( - April_u, - April_bf, - April_ρ, - center_z, +Base.parent(ᶜρ) .= April_ρ +ᶜv = ᶜu +Base.parent(ᶜu) .= April_u +Base.parent(ᶜbf) .= April_bf +April_uforcing = similar(ᶜρ, FT) +April_uforcing .= 0 +April_vforcing = similar(ᶜρ, FT) +April_vforcing .= 0 + +CA.non_orographic_gravity_wave_forcing( + ᶜu, + ᶜv, + ᶜbf, + ᶜρ, + ᶜz, + ᶜlevel, source_level, damp_level, - params.gw_source_ampl, - params.gw_Bw, - params.gw_Bn, - B0, - params.gw_cw, - params.gw_cn, - params.gw_flag, - params.gw_c, - params.gw_c0, - params.gw_nk, + ᶜρ_source, + ᶜu_source, + ᶜv_source, + April_uforcing, + April_vforcing, + gw_ncval, + u_waveforcing, + v_waveforcing, + params, ) +April_uforcing = Base.parent(April_uforcing) fig = generate_empty_figure(); create_plot!( fig; @@ -197,24 +276,36 @@ CairoMakie.save(joinpath(output_dir, "fig6apr.png"), fig) 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 = CA.non_orographic_gravity_wave_forcing( - July_u, - July_bf, - July_ρ, - center_z, +Base.parent(ᶜρ) .= July_ρ +ᶜv = ᶜu +Base.parent(ᶜu) .= July_u +Base.parent(ᶜbf) .= July_bf +July_uforcing = similar(ᶜρ, FT) +July_uforcing .= 0 +July_vforcing = similar(ᶜρ, FT) +July_vforcing .= 0 + +CA.non_orographic_gravity_wave_forcing( + ᶜu, + ᶜv, + ᶜbf, + ᶜρ, + ᶜz, + ᶜlevel, source_level, damp_level, - params.gw_source_ampl, - params.gw_Bw, - params.gw_Bn, - B0, - params.gw_cw, - params.gw_cn, - params.gw_flag, - params.gw_c, - params.gw_c0, - params.gw_nk, + ᶜρ_source, + ᶜu_source, + ᶜv_source, + July_uforcing, + July_vforcing, + gw_ncval, + u_waveforcing, + v_waveforcing, + params, ) +July_uforcing = Base.parent(July_uforcing) + fig = generate_empty_figure(); create_plot!( fig; @@ -228,24 +319,36 @@ CairoMakie.save(joinpath(output_dir, "fig6jul.png"), fig) 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 = CA.non_orographic_gravity_wave_forcing( - Oct_u, - Oct_bf, - Oct_ρ, - center_z, +Base.parent(ᶜρ) .= Oct_ρ +ᶜv = ᶜu +Base.parent(ᶜu) .= Oct_u +Base.parent(ᶜbf) .= Oct_bf +Oct_uforcing = similar(ᶜρ, FT) +Oct_uforcing .= 0 +Oct_vforcing = similar(ᶜρ, FT) +Oct_vforcing .= 0 + +CA.non_orographic_gravity_wave_forcing( + ᶜu, + ᶜv, + ᶜbf, + ᶜρ, + ᶜz, + ᶜlevel, source_level, damp_level, - params.gw_source_ampl, - params.gw_Bw, - params.gw_Bn, - B0, - params.gw_cw, - params.gw_cn, - params.gw_flag, - params.gw_c, - params.gw_c0, - params.gw_nk, + ᶜρ_source, + ᶜu_source, + ᶜv_source, + Oct_uforcing, + Oct_vforcing, + gw_ncval, + u_waveforcing, + v_waveforcing, + params, ) + +Oct_uforcing = Base.parent(Oct_uforcing) fig = generate_empty_figure(); create_plot!( fig;