ConvectiveHeating implemented

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