ConvectiveHeating implemented

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Unreleased
 
+- ConvectiveHeating implemented [#639](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/639)
 - Number format flexibility with set! [#634](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/634)
 - Forcing/drag for primitive models [#635](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/635)
 - RingGrids indexing with leading Colon should now always return another RingGrid instance [#637](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/637)

src/physics/convection.jl

@@ -430,4 +430,68 @@ 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
+
+"""Convective heating as defined by Lee and Kim, 2003, JAS
+implemented as convection parameterization. Fields are
+$(TYPEDFIELDS)"""
+@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...)
+
+# precompute latitudinal mask
+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

test/convection.jl

@@ -4,7 +4,8 @@
 
     for Convection in ( NoConvection,
                         SimplifiedBettsMiller,
-                        DryBettsMiller)
+                        DryBettsMiller,
+                        ConvectiveHeating)
         for Model in (  PrimitiveDryModel,
                         PrimitiveWetModel)