Skip to content

Fix subsidence forcing #15

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mmr0 wants to merge 6 commits into
base: glw/bomex
Choose a base branch
from
Open

Fix subsidence forcing #15

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
4814219
Revert "Revert "- version highlighting the missing step in the subsid…
mmr0 Jun 30, 2025
ae46dfa
Delete random file notes.md
mmr0 Jun 30, 2025
45960cf
attempt at adding plane-average field to forcing using callback
mmr0 Jul 1, 2025
925b7f4
Merge branch 'glw/bomex' of https://github.com/mmr0/Breeze.jl into gl…
mmr0 Jul 1, 2025
df6b138
Update units and formatting in notes.md
navidcy Sep 6, 2025
9af5ee0
Update title for MABL modeling notes
navidcy Sep 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 36 additions & 0 deletions docs/src/notes.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
# Notes on Marine Atmospheric Boundery Layer modeling

## Surface fluxes

### Momentum flux

| Symbol | Name | Unit |
| -------- | ------- |-------- |
| $C^D = 10^{-3}$ | drag coefficient | |
| $G = 0.01$ | gustiness friction velocity | m / s |
| $u^\star = \sqrt(Cᴰ (Δu^2 + Δv^2)) + G$ | friction velocity | m / s |
| $u^{\star 2}u $ | $x$ momentum flux | m³ / s³ |
| $u^{\star 2}v $ | $y$ momentum flux | m³ / s³ |


### Sensible heat flux

| Symbol | Name | Unit |
| -------- | ------- |-------- |
| $C^H = 10^{-3}$ | heat transfer coefficient | |
| $c_p = 1005$ | dry air heat capacity | J / (kg K) |
| $\theta\star = C^H/\sqrt{C^D}\Delta\theta$ | friction temperature | K |
| $J\theta = -u^\star\theta^\star$ | temperature flux | K m / s |
| $J^{SH} = \rho_0 c_p J\theta$ | heat flux | W / m² |


### Latent heat flux

| Symbol | Name | Unit |
| -------- | ------- |-------- |
| $C^v = 10^{-3}$ | vapor transfer coefficient | |
| $L^v = 2.5008 \times 10 ^6$ | latent heat of vaporisation | J / kg |
| $q^\star = Cᵛ / \sqrt{Cᴰ} Δq $ | friction humidity (?) | kg / kg |
| $Jq = -u^\star q^\star$ | specific humidity flux | m / s |
| $J^{LH} = \rho_0 L^v Jq$ | latent heat flux | W / m² |

191 changes: 118 additions & 73 deletions examples/bomex.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ Fq_precip = Forcing(precipitation, field_dependencies=:q, parameters=microphysic

# See Siebesma et al (2003), appendix B1-B2
θ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(8e-3))
q_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(5.2e-5))
q_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1.4e-4))

u★ = 0.28 # m/s
@inline u_drag(x, y, t, u, v, u★) = - u★^2 * u / sqrt(u^2 + v^2)
Expand All @@ -58,12 +58,8 @@ u_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(u_drag, field_depen
v_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(v_drag, field_dependencies=(:u, :v), parameters=u★))

coriolis = FPlane(f=3.76e-5)
uᵍ = Field{Nothing, Nothing, Center}(grid)
vᵍ = Field{Nothing, Nothing, Center}(grid)
uᵍ_bomex = AtmosphericProfilesLibrary.Bomex_geostrophic_u(FT)
vᵍ_bomex = AtmosphericProfilesLibrary.Bomex_geostrophic_v(FT)
set!(uᵍ, z -> uᵍ_bomex(z))
set!(vᵍ, z -> vᵍ_bomex(z))

Fu = Field{Nothing, Nothing, Center}(grid)
Fv = Field{Nothing, Nothing, Center}(grid)
Expand Down Expand Up @@ -99,11 +95,17 @@ end
return - w_dz_Q
end

# These empty fields are being used in the forcing. Missing a step where we calculate the plane averages of the model fields
u_avg = Field{Nothing, Nothing, Center}(grid)
v_avg = Field{Nothing, Nothing, Center}(grid)
θ_avg = Field{Nothing, Nothing, Center}(grid)
q_avg = Field{Nothing, Nothing, Center}(grid)

# ## Test only! see what happens if we use the initial condition as the plane-average profile @ all timesteps
# set!(u_avg, z -> u_bomex(z))
# set!(θ_avg, z -> θ_bomex(z))
# set!(q_avg, z -> q_bomex(z))

wˢ = Field{Nothing, Nothing, Face}(grid)
w_bomex = AtmosphericProfilesLibrary.Bomex_subsidence(FT)
set!(wˢ, z -> w_bomex(z))
Expand All @@ -126,30 +128,32 @@ dqdt_bomex = AtmosphericProfilesLibrary.Bomex_dqtdt(FT)
set!(drying, z -> dqdt_bomex(z))
Fq_drying = Forcing(drying)
q_forcing = (Fq_precip, Fq_drying, Fq_subsidence)
#q_forcing = (Fq_drying, Fq_subsidence)
# q_forcing = Fq_subsidence


Fθ_field = Field{Nothing, Nothing, Center}(grid)
dTdt_bomex = AtmosphericProfilesLibrary.Bomex_dTdt(FT)
set!(Fθ_field, z -> dTdt_bomex(1, z))
Fθ_radiation = Forcing(Fθ_field)
θ_forcing = (Fθ_radiation, Fθ_subsidence)


advection = WENO() #(momentum=WENO(), θ=WENO(), q=WENO(bounds=(0, 1)))
model = NonhydrostaticModel(; grid, advection, buoyancy, coriolis,
tracers = (:θ, :q),
# tracers = (:θ, :q, :qˡ, :qⁱ, :qʳ, :qˢ),
forcing = (q=q_forcing, u=u_forcing, v=v_forcing, θ=θ_forcing),
boundary_conditions = (θ=θ_bcs, q=q_bcs, u=u_bcs))



θϵ = 20
qϵ = 1e-2
θᵢ(x, y, z) = θ_bomex(z) + 1e-2 * θϵ * rand()
qᵢ(x, y, z) = q_bomex(z) + 1e-2 * qϵ * rand()
uᵢ(x, y, z) = u_bomex(z)
set!(model, θ=θᵢ, q=qᵢ, u=uᵢ)

simulation = Simulation(model, Δt=10, stop_time=6hours)
simulation = Simulation(model, Δt=10, stop_time=1hours)
conjure_time_step_wizard!(simulation, cfl=0.7)

T = AquaSkyLES.TemperatureField(model)
Expand Down Expand Up @@ -190,79 +194,120 @@ end

add_callback!(simulation, progress, IterationInterval(10))

u_avg_op = Field(Average(model.velocities.u, dims =(1,2)))

function compute_averages!(model)
u_avg .= u_avg_op
return
end

simulation.callbacks[:update_averages] = Callback(compute_averages!, IterationInterval(1))

# # have a look at the forcing that we're applying
# using Oceananigans.Models: ForcingOperation
# using Oceananigans.Models: ForcingField
# θ_forcing_op = ForcingOperation(:θ, model)
# u_forcing_op = ForcingOperation(:u, model)
# θ_forcing_field = ForcingField(:θ, model)
# u_forcing_field = ForcingField(:u, model)
# compute!(u_forcing_field)
# compute!(θ_forcing_field)

# using Oceananigans.Models: ForcingOperation
# Sʳ = ForcingOperation(:q, model)
# outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★, Sʳ))
outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★))

ow = JLD2Writer(model, outputs,
filename = "bomex.jld2",
# Output slices instead of 3d fields
ow_yz = JLD2Writer(model, outputs,
indices = (1,:,:),
filename = "bomex_yz.jld2",
schedule = TimeInterval(1minutes),
overwrite_existing = true)

simulation.output_writers[:jld2] = ow

ow_xz = JLD2Writer(model, outputs,
indices = (:,1,:),
filename = "bomex_xz.jld2",
schedule = TimeInterval(1minutes),
overwrite_existing = true)

# Output some average profiles
u_prof = Field(Average(model.velocities.u, dims=(1, 2)))
v_prof = Field(Average(model.velocities.v, dims=(1, 2)))
θ_prof = Field(Average(model.tracers.θ, dims=(1, 2)))
q_prof = Field(Average(model.tracers.q, dims=(1, 2)))
qˡ_prof = Field(Average(qˡ, dims=(1, 2)))

ow_profiles = JLD2Writer(model, (; u_prof, v_prof, θ_prof, q_prof, qˡ_prof),
filename = "bomex_profiles.jld2",
schedule = TimeInterval(1minutes),
overwrite_existing = true)



push!(simulation.output_writers, ow_xz, ow_yz, ow_profiles)
run!(simulation)

wt = FieldTimeSeries("bomex.jld2", "w")
θt = FieldTimeSeries("bomex.jld2", "θ")
Tt = FieldTimeSeries("bomex.jld2", "T")
qt = FieldTimeSeries("bomex.jld2", "q")
qˡt = FieldTimeSeries("bomex.jld2", "qˡ")
times = qt.times
Nt = length(θt)

using GLMakie, Printf

fig = Figure(size=(1200, 800), fontsize=12)
axθ = Axis(fig[1, 1], xlabel="x (m)", ylabel="z (m)")
axq = Axis(fig[1, 2], xlabel="x (m)", ylabel="z (m)")
axT = Axis(fig[2, 1], xlabel="x (m)", ylabel="z (m)")
axqˡ = Axis(fig[2, 2], xlabel="x (m)", ylabel="z (m)")
axw = Axis(fig[3, 1], xlabel="x (m)", ylabel="z (m)")

Nt = length(θt)
slider = Slider(fig[4, 1:2], range=1:Nt, startvalue=1)

n = slider.value #Observable(length(θt))
wn = @lift view(wt[$n], :, 1, :)
θn = @lift view(θt[$n], :, 1, :)
qn = @lift view(qt[$n], :, 1, :)
Tn = @lift view(Tt[$n], :, 1, :)
qˡn = @lift view(qˡt[$n], :, 1, :)
title = @lift "t = $(prettytime(times[$n]))"


fig[0, :] = Label(fig, title, fontsize=22, tellwidth=false)

Tmin = 299 #minimum(Tt)
Tmax = 301 #maximum(Tt)
wlim = maximum(abs, wt) / 2
qlim = maximum(abs, qt)
qˡlim = maximum(abs, qˡt) / 2 + 1e-6

Tₛ = θ_bomex(0)
Δθ = θ_bomex(Lz) - θ_bomex(0)
hmθ = heatmap!(axθ, θn, colorrange=(Tₛ, Tₛ+Δθ))
hmq = heatmap!(axq, qn, colorrange=(0, qlim), colormap=:magma)
hmT = heatmap!(axT, Tn, colorrange=(Tmin, Tmax))
hmqˡ = heatmap!(axqˡ, qˡn, colorrange=(0, qˡlim), colormap=:magma)
hmw = heatmap!(axw, wn, colorrange=(-wlim, wlim), colormap=:balance)

# Label(fig[0, 1], "θ", tellwidth=false)
# Label(fig[0, 2], "q", tellwidth=false)
# Label(fig[0, 1], "θ", tellwidth=false)
# Label(fig[0, 2], "q", tellwidth=false)

Colorbar(fig[1, 0], hmθ, label = "θ [K]", vertical=true)
Colorbar(fig[1, 3], hmq, label = "q", vertical=true)
Colorbar(fig[2, 0], hmT, label = "T [K]", vertical=true)
Colorbar(fig[2, 3], hmqˡ, label = "qˡ", vertical=true)
Colorbar(fig[3, 0], hmw, label = "w", vertical=true)

fig

record(fig, "bomex.mp4", 1:Nt, framerate=12) do nn
@info "Drawing frame $nn of $Nt..."
n[] = nn
end
# wt = FieldTimeSeries("bomex.jld2", "w")
# θt = FieldTimeSeries("bomex.jld2", "θ")
# Tt = FieldTimeSeries("bomex.jld2", "T")
# qt = FieldTimeSeries("bomex.jld2", "q")
# qˡt = FieldTimeSeries("bomex.jld2", "qˡ")
# times = qt.times
# Nt = length(θt)

# using GLMakie, Printf

# fig = Figure(size=(1200, 800), fontsize=12)
# axθ = Axis(fig[1, 1], xlabel="x (m)", ylabel="z (m)")
# axq = Axis(fig[1, 2], xlabel="x (m)", ylabel="z (m)")
# axT = Axis(fig[2, 1], xlabel="x (m)", ylabel="z (m)")
# axqˡ = Axis(fig[2, 2], xlabel="x (m)", ylabel="z (m)")
# axw = Axis(fig[3, 1], xlabel="x (m)", ylabel="z (m)")

# Nt = length(θt)
# slider = Slider(fig[4, 1:2], range=1:Nt, startvalue=1)

# n = slider.value #Observable(length(θt))
# wn = @lift view(wt[$n], :, 1, :)
# θn = @lift view(θt[$n], :, 1, :)
# qn = @lift view(qt[$n], :, 1, :)
# Tn = @lift view(Tt[$n], :, 1, :)
# qˡn = @lift view(qˡt[$n], :, 1, :)
# title = @lift "t = $(prettytime(times[$n]))"


# fig[0, :] = Label(fig, title, fontsize=22, tellwidth=false)

# Tmin = 299 #minimum(Tt)
# Tmax = 301 #maximum(Tt)
# wlim = maximum(abs, wt) / 2
# qlim = maximum(abs, qt)
# qˡlim = maximum(abs, qˡt) / 2 + 1e-6

# Tₛ = θ_bomex(0)
# Δθ = θ_bomex(Lz) - θ_bomex(0)
# hmθ = heatmap!(axθ, θn, colorrange=(Tₛ, Tₛ+Δθ))
# hmq = heatmap!(axq, qn, colorrange=(0, qlim), colormap=:magma)
# hmT = heatmap!(axT, Tn, colorrange=(Tmin, Tmax))
# hmqˡ = heatmap!(axqˡ, qˡn, colorrange=(0, qˡlim), colormap=:magma)
# hmw = heatmap!(axw, wn, colorrange=(-wlim, wlim), colormap=:balance)

# # Label(fig[0, 1], "θ", tellwidth=false)
# # Label(fig[0, 2], "q", tellwidth=false)
# # Label(fig[0, 1], "θ", tellwidth=false)
# # Label(fig[0, 2], "q", tellwidth=false)

# Colorbar(fig[1, 0], hmθ, label = "θ [K]", vertical=true)
# Colorbar(fig[1, 3], hmq, label = "q", vertical=true)
# Colorbar(fig[2, 0], hmT, label = "T [K]", vertical=true)
# Colorbar(fig[2, 3], hmqˡ, label = "qˡ", vertical=true)
# Colorbar(fig[3, 0], hmw, label = "w", vertical=true)

# fig

# record(fig, "bomex.mp4", 1:Nt, framerate=12) do nn
# @info "Drawing frame $nn of $Nt..."
# n[] = nn
# end