Skip to content

Commit

Permalink
Merge pull request #198 from SpeedyWeather/mk/vertadvection
Browse files Browse the repository at this point in the history 
Vertical advection of u,v + tendencies (vor, pres grads)
  • Loading branch information
@milankl
milankl committed Dec 2, 2022
2 parents afaad59 + 35fdc8e commit 68056ef
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 73 deletions.
15 changes: 10 additions & 5 deletions src/diagnostic_variables.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,10 @@ struct Tendencies{NF<:AbstractFloat,Grid<:AbstractGrid{NF}}
div_tend ::LowerTriangularMatrix{Complex{NF}} # Divergence of horizontal wind field [1/s]
temp_tend ::LowerTriangularMatrix{Complex{NF}} # Absolute temperature [K]
humid_tend::LowerTriangularMatrix{Complex{NF}} # Specific humidity [g/kg]
u_tend ::Grid # zonal velocity
v_tend ::Grid # meridional velocity
u_tend ::LowerTriangularMatrix{Complex{NF}} # zonal velocity (spectral)
v_tend ::LowerTriangularMatrix{Complex{NF}} # meridional velocity (spectral)
u_tend_grid ::Grid # zonal velocity (grid)
v_tend_grid ::Grid # meridinoal velocity (grid)
humid_grid_tend ::Grid # specific humidity
temp_grid_tend ::Grid # temperature
end
Expand All @@ -25,13 +27,16 @@ function Base.zeros(::Type{Tendencies},
div_tend = zeros(LowerTriangularMatrix{Complex{NF}},lmax+2,mmax+1) # divergence
temp_tend = zeros(LowerTriangularMatrix{Complex{NF}},lmax+2,mmax+1) # absolute Temperature
humid_tend = zeros(LowerTriangularMatrix{Complex{NF}},lmax+2,mmax+1) # specific humidity
u_tend = zeros(Grid{NF},nresolution) # zonal velocity
v_tend = zeros(Grid{NF},nresolution) # meridional velocity
u_tend = zeros(LowerTriangularMatrix{Complex{NF}},lmax+2,mmax+1) # zonal velocity
v_tend = zeros(LowerTriangularMatrix{Complex{NF}},lmax+2,mmax+1) # meridional velocity
u_tend_grid = zeros(Grid{NF},nresolution) # zonal velocity
v_tend_grid = zeros(Grid{NF},nresolution) # meridional velocity
humid_grid_tend = zeros(Grid{NF},nresolution) # specific humidity
temp_grid_tend = zeros(Grid{NF},nresolution) # temperature

return Tendencies( vor_tend,div_tend,temp_tend,humid_tend,
u_tend,v_tend,humid_grid_tend,temp_grid_tend)
u_tend,v_tend,u_tend_grid,v_tend_grid,
humid_grid_tend,temp_grid_tend)
end

"""
Expand Down
17 changes: 17 additions & 0 deletions src/tendencies.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -393,6 +393,15 @@ function get_spectral_tendencies!(Prog::PrognosticVariables{NF},

end

"""
add_tendencies!(tend::LowerTriangularMatrix{NF}, # tendency to accumulate into
term1::LowerTriangularMatrix{NF}, # with term1
term2::LowerTriangularMatrix{NF} # and term2
) where NF # number format real or complex
Accumulates three `LowerTriangularMatrix`s element-wise into the first.
tend += term1 + term2."""
function add_tendencies!( tend::LowerTriangularMatrix{NF}, # tendency to accumulate into
term1::LowerTriangularMatrix{NF}, # with term1
term2::LowerTriangularMatrix{NF} # and term2
Expand All @@ -403,6 +412,14 @@ function add_tendencies!( tend::LowerTriangularMatrix{NF}, # tendency to ac
end
end

"""
add_tendencies!(tend::LowerTriangularMatrix{NF}, # tendency to accumulate into
term::LowerTriangularMatrix{NF}, # with term
) where NF # number format real or complex
Accumulates two `LowerTriangularMatrix`s element-wise into the first.
tend += term."""
function add_tendencies!( tend::LowerTriangularMatrix{NF}, # tendency to accumulate into
term::LowerTriangularMatrix{NF} # with term
) where NF # number format real or complex
Expand Down
140 changes: 72 additions & 68 deletions src/tendencies_dynamics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ function vertical_velocity!(Diag::DiagnosticVariables{NF},
V = Diag.layers[k].grid_variables.V_grid
D = Diag.layers[k].grid_variables.div_grid

# σ_tend/σ_m sits on half layers below (k+1/2), but its 0 at
# σ_tend & σ_m sit on half layers below (k+1/2), but its 0 at
# k=1/2 and nlev+1/2, don't explicitly store k=1/2
σ_tend = Diag.layers[k].dynamics_variables.σ_tend # actually on half levels
σ_m = Diag.layers[k].dynamics_variables.σ_tend # actually on half levels
Expand All @@ -112,6 +112,7 @@ function vertical_velocity!(Diag::DiagnosticVariables{NF},
σ_tend_below = Diag.layers[kmax].dynamics_variables.σ_tend
σ_m_below = Diag.layers[kmax].dynamics_variables.σ_m

# TODO check whether coslat unscaling is needed
for ij in eachgridpoint(U,V,D,U_mean,V_mean,div_mean,dpres_dlon_grid,dpres_dlat_grid)
uv∇p_ij = (U[ij]-U_mean[ij])*dpres_dlon_grid[ij] + (V[ij]-V_mean[ij])*dpres_dlat_grid[ij]
uv∇p[ij] = uv∇p_ij
Expand All @@ -129,6 +130,47 @@ function vertical_velocity!(Diag::DiagnosticVariables{NF},
end
end

function vertical_advection!( diagn::DiagnosticVariables,
model::PrimitiveEquationModel)

@unpack σ_levels_thick⁻¹_half, nlev = model.geometry
@boundscheck nlev == diagn.nlev || throw(BoundsError)

# ALL LAYERS (but use indexing tricks to avoid out of bounds access for top/bottom)
@inbounds for k in 1:nlev
# for k==1 "above" term is 0, for k==nlev "below" term is zero
# avoid out-of-bounds indexing with k_above, k_below as follows
k_above = k == 1 ? nlev : k-1 # wrap around to access M_nlev+1/2 = 0 (which zeros that term)
k_below = min(k+1,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term)

# mass fluxes, M_1/2 = M_nlev+1/2 = 0, but k=1/2 isn't explicitly stored
σ_tend_above = diagn.layers[k_above].dynamics_variables.σ_tend
σ_tend_below = diagn.layers[k].dynamics_variables.σ_tend

# zonal wind
u_tend = diagn.layers[k].tendencies.u_tend_grid
U_above = diagn.layers[k_above].grid_variables.U_grid
U = diagn.layers[k].grid_variables.U_grid
U_below = diagn.layers[k_below].grid_variables.U_grid

# meridional wind
v_tend = diagn.layers[k].tendencies.v_tend_grid
V_above = diagn.layers[k_above].grid_variables.V_grid
V = diagn.layers[k].grid_variables.V_grid
V_below = diagn.layers[k_below].grid_variables.V_grid

Δσk = σ_levels_thick⁻¹_half[k] # = 1/(2Δp_k), for convenience

# TODO check whether coslat unscaling is needed
@inbounds for ij in eachgridpoint(u_tend,v_tend)
u_tend[ij] = (σ_tend_above[ij]*(U_above[ij] - U[ij]) +
σ_tend_below[ij]*(U[ij] - U_below[ij]))*Δσk
v_tend[ij] = (σ_tend_above[ij]*(V_above[ij] - V[ij]) +
σ_tend_below[ij]*(V[ij] - V_below[ij]))*Δσk
end
end
end

"""
Compute the temperature anomaly in grid point space
"""
Expand All @@ -148,81 +190,43 @@ function temperature_grid_anomaly!(Diag::DiagnosticVariables{NF}, # Diagnostic v

end

"""
Compute the spectral tendency of the zonal wind
"""
function zonal_wind_tendency!(Diag::DiagnosticVariables{NF}, # Diagnostic variables
M
)where {NF<:AbstractFloat}

@unpack u_tend = Diag.tendencies
@unpack u_grid,v_grid,vor_grid,temp_grid_anomaly= Diag.grid_variables
@unpack sigma_tend,pres_surf_gradient_grid_x,pres_surf_gradient_grid_y,sigma_u = Diag.intermediate_variables


@unpack rgas,σ_levels_half⁻¹_2 = M.GeoSpectral.geometry #I think this is dhsr

_,_,nlev = size(u_grid)


#Update px,py

pres_surf_gradient_grid_x = rgas*pres_surf_gradient_grid_x
pres_surf_gradient_grid_y = rgas*pres_surf_gradient_grid_y



for k in 2:nlev
sigma_u[:,:,k] = sigma_tend[:,:,k].*(u_grid[:,:,k] - u_grid[:,:,k-1])
end


for k in 1:nlev
u_tend[:,:,k] = u_tend[:,:,k] + v_grid[:,:,k].*vor_grid[:,:,k]
- temp_grid_anomaly[:,:,k].*pres_surf_gradient_grid_x
- (sigma_u[:,:,k+1] + sigma_u[:,:,k])*σ_levels_half⁻¹_2[k]
end




end



"""
Compute the spectral tendency of the meridional wind
"""
function meridional_wind_tendency!(Diag::DiagnosticVariables{NF}, # Diagnostic variables
M
)where {NF<:AbstractFloat}

@unpack v_tend = Diag.tendencies
@unpack vor_grid,u_grid,v_grid,temp_grid_anomaly =Diag.grid_variables
@unpack sigma_tend,sigma_u,pres_surf_gradient_grid_x,pres_surf_gradient_grid_y = Diag.intermediate_variables

function uv_tendencies!(diagn::DiagnosticVariablesLayer,
surf::SurfaceVariables,
model::PrimitiveEquationModel)

@unpack rgas,σ_levels_half⁻¹_2 = M.GeoSpectral.geometry #I think this is dhsr
@unpack f_coriolis, coslat⁻² = model.geometry
@unpack R_dry = model.constants

@unpack u_tend_grid, v_tend_grid = diagn.tendencies # already contains vertical advection
U = diagn.grid_variables.U_grid # U = u*coslat
V = diagn.grid_variables.V_grid # V = v*coslat
vor = diagn.grid_variables.vor_grid # relative vorticity
dpres_dx = surf.dpres_dlon_grid # zonal gradient of logarithm of surface pressure
dpres_dy = surf.dpres_dlat_grid # meridional gradient thereof
Tᵥ = diagn.grid_variables.temp_virt_grid # virtual temperature

_,_,nlev = size(u_grid)
# precompute ring indices and boundscheck
rings = eachring(u_tend_grid,v_tend_grid,U,V,vor,dpres_dx,dpres_dy,Tᵥ)


for k in 2:nlev
sigma_u[:,:,k] = sigma_tend[:,:,k].*(v_grid[:,:,k] - v_grid[:,:,k-1])
end



for k in 1:nlev
v_tend[:,:,k] = v_tend[:,:,k] + u_grid[:,:,k].*vor_grid[:,:,k]
- temp_grid_anomaly[:,:,k].*pres_surf_gradient_grid_y
- (sigma_u[:,:,k+1] + sigma_u[:,:,k])*σ_levels_half⁻¹_2[k]
@inbounds for (j,ring) in enumerate(rings)
coslat⁻²j = coslat⁻²[j]
f = f_coriolis[j]
for ij in ring
ω = vor[ij] + f # absolute vorticity
RTᵥ = R_dry*Tᵥ[ij] # gas constant (dry air) times virtual temperature
u_tend_grid[ij] = (u_tend_grid[ij] + V[ij]*ω - RTᵥ*dpres_dx[ij])*coslat⁻²j
v_tend_grid[ij] = (v_tend_grid[ij] - U[ij]*ω - RTᵥ*dpres_dy[ij])*coslat⁻²j
end
end

# convert to spectral space
@unpack u_tend, v_tend = diagn.tendencies # spectral fields
S = model.spectral_transform

spectral!(u_tend,u_tend_grid,S)
spectral!(v_tend,v_tend_grid,S)


return nothing
end

"""
Expand Down

0 comments on commit 68056ef

Please sign in to comment.