Skip to content

Commit

Permalink
revise 1
Browse files Browse the repository at this point in the history
  • Loading branch information
谢萧涯 authored and 谢萧涯 committed Jul 12, 2024
1 parent e2bf71b commit 5750311
Show file tree
Hide file tree
Showing 5 changed files with 80 additions and 50 deletions.
2 changes: 2 additions & 0 deletions artifacts/Artifacts.toml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
[topo-info]
git-tree-sha1 = "9970989c13c2ad015ab5738d42fd9f5b320f1451"
122 changes: 75 additions & 47 deletions src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,19 +28,19 @@ 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})
ᶜlevel = similar(Y.c.ρ,FT)
for i in 1: size(Fields.field_values(Y.c.ρ),4)
fill!(Fields.level(ᶜlevel,i),i);
source_level_z = similar(Fields.level(Y.c.ρ,1), Tuple{FT,FT})
ᶜlevel = similar(Y.c.ρ, FT)
for i in 1:size(Fields.field_values(Y.c.ρ), 4)
fill!(Fields.level(ᶜlevel,i), i)
end
damp_level_z=source_level_z
damp_level_z= source_level_z
return (;
gw_source_height = source_height,
gw_source_ampl = Bt_0 .* ones(FT, axes(Fields.level(Y.c.ρ, 1))),
gw_Bw = Bw .* ones(FT, axes(Fields.level(Y.c.ρ, 1))),
gw_Bn = Bn .* ones(FT, axes(Fields.level(Y.c.ρ, 1))),
gw_B0 = similar(c),
gw_c = c,
gw_c = c,
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,
Expand All @@ -53,10 +53,10 @@ function non_orographic_gravity_wave_cache(
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.ρ),
uforcing=similar(Y.c.ρ),
vforcing=similar(Y.c.ρ),
u_phy = similar(Y.c.ρ),
v_phy = similar(Y.c.ρ),
uforcing = similar(Y.c.ρ),
vforcing = similar(Y.c.ρ),
)
end

Expand All @@ -82,6 +82,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})
ᶜlevel = similar(Y.c.ρ, FT)
for i in 1:size(Fields.field_values(Y.c.ρ), 4)
fill!(Fields.level(ᶜlevel,i), i);
end
damp_level_z = source_level_z

# This is GFDL source specs -> a smooth function
# source_ampl = @. Bt_0 +
# Bt_n * FT(0.5) * (FT(1) + tanh((lat - ϕ0_n) / dϕ_n)) +
Expand Down Expand Up @@ -123,10 +130,11 @@ function non_orographic_gravity_wave_cache(
damp_level_z,
source_level = similar(Fields.level(Y.c.ρ, 1)),
damp_level = similar(Fields.level(Y.c.ρ, 1)),
u_phy=similar(Y.c.ρ),
v_phy=similar(Y.c.ρ),
uforcing=similar(Y.c.ρ),
vforcing=similar(Y.c.ρ),
ᶜlevel,
u_phy = similar(Y.c.ρ),
v_phy = similar(Y.c.ρ),
uforcing = similar(Y.c.ρ),
vforcing = similar(Y.c.ρ),
)
end

Expand All @@ -141,7 +149,19 @@ function non_orographic_gravity_wave_tendency!(
(; ᶜT,) = p.core
(; ᶜts) = p.precomputed
(; params) = p
(; ᶜdTdz, ᶜbuoyancy_frequency,source_level,source_level_z,damp_level,source_level_z,u_phy,v_phy,uforcing,vforcing,ᶜlevel) = p.non_orographic_gravity_wave
(;
ᶜdTdz,
ᶜbuoyancy_frequency,
source_level,
source_level_z,
damp_level,
damp_level_z,
u_phy,
v_phy,
uforcing,
vforcing,
ᶜlevel,
) = p.non_orographic_gravity_wave
(; model_config) = p.atmos
(;
gw_source_ampl,
Expand Down Expand Up @@ -184,29 +204,44 @@ function non_orographic_gravity_wave_tendency!(
if model_config isa SingleColumnModel
# source level: the index of the level that is closest to the source height

Operators.column_mapreduce!(reduce_fun1, source_level_z, ᶜz, ᶜlevel) do z, level
(abs.(z .- gw_source_height),level)
Operators.column_mapreduce!(
reduce_fun1,
source_level_z,
ᶜz,
ᶜlevel,
) do z, level
(abs.(z .- gw_source_height), level)
end
source_level=source_level_z.:2
source_level = source_level_z.:2

Operators.column_mapreduce!(sign, +, damp_level, ᶜz)

elseif model_config isa SphericalModel
(; ᶜp) = p.precomputed
# source level: the index of the highest level whose pressure is higher than source pressure

Operators.column_mapreduce!(reduce_fun2, source_level_z, ᶜp, ᶜlevel) do p, level
(p .- gw_source_pressure,level)
Operators.column_mapreduce!(
reduce_fun2,
source_level_z,
ᶜp,
ᶜlevel
) do p, level
(p .- gw_source_pressure, level)
end
source_level=source_level_z.:2
source_level = source_level_z.:2


# damp level: the index of the lowest level whose pressure is lower than the damp pressure

Operators.column_mapreduce!(reduce_fun3, damp_level_z, ᶜp, ᶜlevel) do p, level
(p .- gw_damp_pressure,level)
Operators.column_mapreduce!(
reduce_fun3,
damp_level_z,
ᶜp,
ᶜlevel
) do p, level
(p .- gw_damp_pressure, level)
end
damp_level=damp_level_z.:2
damp_level = damp_level_z.:2

end

Expand Down Expand Up @@ -287,17 +322,20 @@ function non_orographic_gravity_wave_forcing(
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])
ᶜ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

B0 =@. sign(c_hat0)* (Bw * exp(-log(2.0) * ((c * flag + c_hat0 * (1 - flag) - c0) / cw)^2)+ Bn *exp(-log(2.0) * ((c * flag + c_hat0 * (1 - flag) - c0) / cn)^2))
B0 = @. sign(c_hat0) * (
Bw * exp(-log(2.0) * ((c * flag + c_hat0 * (1 - flag) - c0) / cw)^2) +
Bn * exp(-log(2.0) * ((c * flag + c_hat0 * (1 - flag) - c0) / cn)^2)
)

Bsum = sum(abs.(B0))
if (Bsum == 0.0)
Expand Down Expand Up @@ -415,30 +453,20 @@ function calc_intermitency(ρ_source_level, source_ampl, nk, Bsum)
return (source_ampl / ρ_source_level / nk) / Bsum
end

function reduce_fun1(a,b)
return ifelse(
a[1]<b[1],
a,
b,
)

function reduce_fun1(a, b)
return ifelse(a[1] < b[1], a, b)
end

function reduce_fun2(a,b)
return ifelse(
b[1]<0,
a,
b,
)
function reduce_fun2(a, b)
return ifelse(b[1] < 0, a, b)
end

function reduce_fun3(a,b)
return ifelse(
a[1]>0,
b,
a,
)
function reduce_fun3(a, b)
return ifelse(a[1] > 0, b, a)
end





2 changes: 1 addition & 1 deletion test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,7 @@ 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)
B0 = similar(params.gw_c)

# Jan
month = Dates.month.(time)
Expand Down
2 changes: 1 addition & 1 deletion test/parameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_mima.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ 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)
B0 = similar(params.gw_c)
# nogw forcing

kmax = length(pfull) - 1
Expand Down
2 changes: 1 addition & 1 deletion ...ameterized_tendencies/gravity_wave/non_orographic_gravity_wave/nogw_test_single_column.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ ENV["GKSwstype"] = "nul"
output_dir = "nonorographic_gravity_wave_test_single_column"
mkpath(output_dir)

B0=similar(params.gw_c)
B0 = similar(params.gw_c)

# Jan
Jan_u = mean(center_u_mean[:, month .== 1], dims = 2)[:, 1]
Expand Down

0 comments on commit 5750311

Please sign in to comment.