diff --git a/src/physics/convection.jl b/src/physics/convection.jl index 44c57421b..a0a0bc6dc 100644 --- a/src/physics/convection.jl +++ b/src/physics/convection.jl @@ -430,4 +430,63 @@ function dry_adiabat!( # level of zero buoyancy is reached when the loop stops, but in case it's at the top it's still buoyant level_zero_buoyancy = k + (1-buoyant) return level_zero_buoyancy +end + +export ConvectiveHeating +@kwdef struct ConvectiveHeating{NF} <: AbstractConvection + # DIMENSION + nlat::Int + + "[OPTION] Q_max heating strength as 1K/time_scale" + time_scale::Second = Hour(12) + + "[OPTION] Pressure of maximum heating [hPa]" + p₀::NF = 525 + + "[OPTION] Vertical extent of heating [hPa]" + σₚ::NF = 200 + + "[OPTION] Latitude of heating [˚N]" + θ₀::NF = 0 + + "[OPTION] Latitudinal width of heating [˚]" + σθ::NF = 20 + + # precomputed latitude mask + lat_mask::Vector{NF} = zeros(NF, nlat) +end + +# generator +ConvectiveHeating(SG::SpectralGrid; kwargs...) = ConvectiveHeating{SG.NF}(nlat=SG.nlat; kwargs...) + +function initialize!(C::ConvectiveHeating, model::PrimitiveEquation) + + (; latd) = model.geometry + (; θ₀, σθ) = C + + for (j, θ) in enumerate(latd) + # Lee and Kim, 2003, eq. 2 + C.lat_mask[j] = cosd((θ-θ₀)/σθ)^2 + end +end + +function convection!( + column::ColumnVariables, + scheme::ConvectiveHeating, + model::PrimitiveEquation, +) + # escape immediately if not in the tropics + abs(column.latd) > scheme.σθ && return nothing + + p₀ = scheme.p₀*100 # hPa -> Pa + σₚ = scheme.σₚ*100 # hPa -> Pa + cos²θ_term = scheme.lat_mask[column.jring] + Qmax = 1/Second(scheme.time_scale).value + + for k in eachindex(column) + p = column.pres[k] # Pressure in Pa + + # Lee and Kim, 2003, eq. 2, + column.temp_tend[k] += Qmax*exp(-((p-p₀)/σₚ)^2 / 2)*cos²θ_term + end end \ No newline at end of file