Skip to content

Commit

Permalink
foo6
Browse files Browse the repository at this point in the history
  • Loading branch information
谢萧涯 authored and 谢萧涯 committed Sep 10, 2024
1 parent 1b4266b commit 15a7648
Showing 1 changed file with 17 additions and 19 deletions.
36 changes: 17 additions & 19 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 @@ -314,7 +314,7 @@ function non_orographic_gravity_wave_forcing(
v_waveforcing,
p,
) where {nc}
#unpack parameters
# unpack parameters
(;
gw_source_ampl,
gw_Bw,
Expand All @@ -336,7 +336,7 @@ function non_orographic_gravity_wave_forcing(

FT = eltype(ᶜρ) # Define the floating point type

#generate the field of ρ,u,v,z with on level shifted up
# Using interpolate operator, generate the field of ρ,u,v,z with on level shifted up
ρ_endlevel = Fields.level(ᶜρ, Spaces.nlevels(axes(ᶜρ)))
ρ_endlevel_m1 = Fields.level(ᶜρ, Spaces.nlevels(axes(ᶜρ)) - 1)
Boundary_value = Fields.Field(
Expand Down Expand Up @@ -381,7 +381,7 @@ function non_orographic_gravity_wave_forcing(
#We use StaticBitVector here because the unrolled_reduce function in Julia can cause memory allocation issues when the mask has more than 32 elements。
#StaticBitVector stores 8 boolean values in a UInt8, allowing efficient storage for up to 256 gravity wave break data.
level_end = Spaces.nlevels(axes(ᶜρ))

# loop over all wave lengths
for ink in 1:gw_nk

Expand Down Expand Up @@ -427,9 +427,9 @@ function non_orographic_gravity_wave_forcing(
)

# Accumulate zonal wave forcing in every column
waveforcing_column_accumulate!(u_waveforcing,mask_u,input_u,gw_c,ink,gw_nk,level_end,gw_c0,gw_ncval)
waveforcing_column_accumulate!(u_waveforcing,mask_u,input_u,gw_c,gw_c0,gw_nk,ink,level_end,gw_ncval)
# Accumulate meridional wave forcing in every column
waveforcing_column_accumulate!(v_waveforcing,mask_v,input_v,gw_c,ink,gw_nk,level_end,gw_c0,gw_ncval)
waveforcing_column_accumulate!(v_waveforcing,mask_v,input_v,gw_c,gw_c0,gw_nk,ink,level_end,gw_ncval)

#extract the momentum flux outside the model top.
u_waveforcing_top = p.scratch.temp_field_level
Expand Down Expand Up @@ -462,11 +462,11 @@ function non_orographic_gravity_wave_forcing(
0,
)

#interpolate the waveforcing from center to face
# interpolate the waveforcing from center to face
gw_average!(u_waveforcing, p.scratch.ᶜtemp_scalar)
gw_average!(v_waveforcing, p.scratch.ᶜtemp_scalar)

#The momentum flux outside the model top will be evenly deposited onto the levels between the damp level and the model top.
# The momentum flux outside the model top will be evenly deposited onto the levels between the damp level and the model top.
@. u_waveforcing = gw_deposit(
u_waveforcing_top,
u_waveforcing,
Expand All @@ -482,17 +482,19 @@ function non_orographic_gravity_wave_forcing(
level_end,
)

#update gravity wave forcing
# update gravity wave forcing
@. uforcing = uforcing + u_waveforcing
@. vforcing = vforcing + v_waveforcing

end
return nothing
end

# calculate the intermittency factor eps -> assuming constant Δc.
function waveforcing_column_accumulate!(waveforcing,mask,input,c,ink,nk,level_end,c0,gw_ncval::Val{nc}) where {nc}
# Using column_accumulate function, calculate the gravity wave forcing at each point.
function waveforcing_column_accumulate!(waveforcing,mask,input,c,c0,nk,ink,level_end,gw_ncval::Val{nc}) where {nc}
FT=eltype(waveforcing)
# Here we use column_accumulate function to pass the variable B0 and mask through different levels, and calculate waveforcing
# at each level.
Operators.column_accumulate!(
waveforcing,
input;
Expand All @@ -519,14 +521,15 @@ function waveforcing_column_accumulate!(waveforcing,mask,input,c,ink,nk,level_en
)

FT1 = typeof(u_kp1)
kwv = 2.0 * π / ((30.0 * (10.0^ink)) * 1.e3) #wave number of gravity waves
kwv = 2.0 * π / ((30.0 * (10.0^ink)) * 1.e3) # wave number of gravity waves
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)

# calculate momentum flux carried by gravity waves with different phase speeds.
B0, Bsum = if level == 1
mask = StaticBitVector{nc}(_ -> true)
B1 = ntuple(
Expand Down Expand Up @@ -560,14 +563,12 @@ function waveforcing_column_accumulate!(waveforcing,mask,input,c,ink,nk,level_en
else
B0_or_NaNs, Bsum_or_NaN
end
#calculate momentum flux carried by gravity waves with different phase speeds.

if level >= source_level - 1
fm = FT1(0.0)
#We use the unrolled_reduce function here because it performs better for parallel execution on the GPU. Additionally, it also helps speed up the setindex function
# We use the unrolled_reduce function here because it performs better for parallel execution on the GPU, avoiding type instabilities.
(mask, fm) = unrolled_reduce(
Val(nc),
(mask, fm),
(mask, FT1(0.0)),
) do (mask, fm), (n)
if (mask[n]) == true
c_hat = c[n] - u_kp1 # c0mu
Expand Down Expand Up @@ -633,6 +634,7 @@ function waveforcing_column_accumulate!(waveforcing,mask,input,c,ink,nk,level_en
end
end

# calculate the intermittency factor eps -> assuming constant Δc.
function calc_intermitency(ρ_source_level, source_ampl, nk, Bsum)
return (source_ampl / ρ_source_level / nk) / Bsum
end
Expand All @@ -653,10 +655,6 @@ function gw_deposit(wave_forcing_top, wave_forcing, damp_level, level, height)
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(;)
Expand Down

0 comments on commit 15a7648

Please sign in to comment.