Merge pull request #439 from SpeedyWeather/mk/condensation

ImmediateCondensation and Clausius-Clapeyron

src/abstract_types.jl

@@ -30,20 +30,11 @@ abstract type AbstractColumnVariables{NF} end
 abstract type AbstractForcing{NF} end
 abstract type AbstractDrag{NF} end
 
-# VERTICAL ADVECTION (PrimitiveEquation)
-abstract type VerticalAdvection{NF,B} end
-
-# SOLAR RADIATION
-abstract type AbstractSolarDeclination{NF} end
-abstract type AbstractSolarTimeCorrection{NF} end
-abstract type AbstractZenith{NF,Grid} end
-
 # PARAMETERIZATIONS
 abstract type AbstractParameterization{NF} end
 abstract type BoundaryLayerDrag{NF} <: AbstractParameterization{NF} end
 abstract type TemperatureRelaxation{NF} <: AbstractParameterization{NF} end
 abstract type VerticalDiffusion{NF} <: AbstractParameterization{NF} end
-abstract type AbstractThermodynamics{NF} <: AbstractParameterization{NF} end
 abstract type AbstractCondensation{NF} <: AbstractParameterization{NF} end
 abstract type AbstractSurfaceWind{NF} <: AbstractParameterization{NF} end
 abstract type AbstractSurfaceThermodynamics{NF} <: AbstractParameterization{NF} end

src/dynamics/initial_conditions.jl

@@ -455,7 +455,7 @@ function initialize_humidity!(  progn::PrognosticVariables,
     temp_grid = gridded(progn.layers[end].timesteps[1].temp,model.spectral_transform)
     humid_surf_grid = zero(pres_surf_grid)
     for ij in eachgridpoint(humid_surf_grid)
-        q_sat = saturation_humidity(temp_grid[ij],exp(pres_surf_grid[ij]),model.thermodynamics)
+        q_sat = saturation_humidity(temp_grid[ij],exp(pres_surf_grid[ij]),model.clausis_clapeyron)
         humid_surf_grid[ij] = relhumid_ref*q_sat
     end
 

src/dynamics/models.jl

@@ -165,7 +165,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevic
     physics::Bool = true
     boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
     temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid)
-    static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid)
+    static_energy_diffusion::VerticalDiffusion{NF} = NoVerticalDiffusion(spectral_grid)
     surface_thermodynamics::AbstractSurfaceThermodynamics{NF} = SurfaceThermodynamicsConstant(spectral_grid)
     surface_wind::AbstractSurfaceWind{NF} = SurfaceWind(spectral_grid)
     surface_heat_flux::AbstractSurfaceHeat{NF} = SurfaceSensibleHeat(spectral_grid)
@@ -255,12 +255,12 @@ Base.@kwdef mutable struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevic
 
     # PHYSICS/PARAMETERIZATIONS
     physics::Bool = true
-    thermodynamics::Thermodynamics{NF} = Thermodynamics(spectral_grid,atmosphere)
+    clausis_clapeyron::AbstractClausiusClapeyron{NF} = ClausiusClapeyron(spectral_grid,atmosphere)
     boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
     temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid)
-    static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid)
-    humidity_diffusion::VerticalDiffusion{NF} = HumidityDiffusion(spectral_grid)
-    large_scale_condensation::AbstractCondensation{NF} = SpeedyCondensation(spectral_grid)
+    static_energy_diffusion::VerticalDiffusion{NF} = NoVerticalDiffusion(spectral_grid)
+    humidity_diffusion::VerticalDiffusion{NF} = NoVerticalDiffusion(spectral_grid)
+    large_scale_condensation::AbstractCondensation{NF} = ImmediateCondensation(spectral_grid)
     surface_thermodynamics::AbstractSurfaceThermodynamics{NF} = SurfaceThermodynamicsConstant(spectral_grid)
     surface_wind::AbstractSurfaceWind{NF} = SurfaceWind(spectral_grid)
     surface_heat_flux::AbstractSurfaceHeat{NF} = SurfaceSensibleHeat(spectral_grid)

src/dynamics/vertical_advection.jl

@@ -1,3 +1,5 @@
+abstract type VerticalAdvection{NF,B} end
+
 # Dispersive and diffusive advection schemes `NF` is the type, `B` the half-stencil size
 abstract type DiffusiveVerticalAdvection{NF, B}  <: VerticalAdvection{NF, B} end
 abstract type DispersiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end

src/physics/column_variables.jl

@@ -141,4 +141,5 @@ function reset_column!(column::ColumnVariables{NF}) where NF
 end
 
 # iterator for convenience
-eachlayer(column::ColumnVariables) = Base.OneTo(column.nlev)
+eachlayer(column::ColumnVariables) = Base.OneTo(column.nlev)
+Base.eachindex(column::ColumnVariables) = Base.OneTo(column.nlev)

src/physics/define_column.jl

@@ -63,7 +63,6 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     surface_air_density::NF = 0
     const sat_humid::Vector{NF} = zeros(NF,nlev)                # Saturation specific humidity [kg/kg]
     const rel_humid::Vector{NF} = zeros(NF,nlev)                # Relative humidity [1]
-    const sat_vap_pres::Vector{NF} = zeros(NF,nlev)             # Saturation vapour pressure [Pa]
     const dry_static_energy::Vector{NF} = zeros(NF,nlev)        # Dry static energy
     const moist_static_energy::Vector{NF} = zeros(NF,nlev)      # Moist static energy
     const sat_moist_static_energy::Vector{NF} = zeros(NF,nlev)  # Saturation moist static energy

src/physics/large_scale_condensation.jl

@@ -76,12 +76,21 @@ function large_scale_condensation!(
     return nothing
 end
 
-"""Function barrier only."""
+# function barrier for all AbstractCondensation
 function large_scale_condensation!( 
     column::ColumnVariables,
     model::PrimitiveWet,
 )
-    large_scale_condensation!(column,model.large_scale_condensation,
+    large_scale_condensation!(column,model.large_scale_condensation,model)
+end
+
+# function barrier for SpeedyCondensation to unpack model
+function large_scale_condensation!( 
+    column::ColumnVariables,
+    scheme::SpeedyCondensation,
+    model::PrimitiveWet,
+)
+    large_scale_condensation!(column,scheme,
         model.geometry,model.constants,model.time_stepping)
 end
 
@@ -143,6 +152,94 @@ function large_scale_condensation!(
             ΔpₖΔt_gρ = σ_levels_thick[k]*pₛΔt_gρ                    # Formula 4 *Δt for [m] of rain during Δt
             column.precip_large_scale += -ΔpₖΔt_gρ * humid_tend_k   # Formula 25, unit [m]
 
+            # only accumulate into humid_tend now to allow humid_tend != 0 before this scheme is called
+            humid_tend[k] += humid_tend_k
+        end
+    end
+end
+
+"""
+Large scale condensation as with immediate precipitation.
+$(TYPEDFIELDS)"""
+Base.@kwdef struct ImmediateCondensation{NF<:AbstractFloat} <: AbstractCondensation{NF}
+    "Flux limiter for latent heat release [K] per timestep"
+    max_heating::NF = 0.2
+
+    "Latent heat of evaporation divided by specific heat of dry air "
+    latent_heat_cₚ::Base.RefValue{NF} = Ref(zero(NF))
+end
+
+ImmediateCondensation(SG::SpectralGrid;kwargs...) = ImmediateCondensation{SG.NF}(;kwargs...)
+
+"""
+$(TYPEDSIGNATURES)
+Initialize the SpeedyCondensation scheme."""
+function initialize!(scheme::ImmediateCondensation,model::PrimitiveEquation)
+    # for latent heat release dT/dt = -(Lᵥ/cₚ)*dq/dt
+    scheme.latent_heat_cₚ[] = model.atmosphere.latent_heat_condensation/    # [J/kg]
+                                model.atmosphere.cₚ     # [J/K/kg] specific heat capacity, const pres
+end
+
+# function barrier for ImmediateCondensation to unpack model
+function large_scale_condensation!( 
+    column::ColumnVariables,
+    scheme::ImmediateCondensation,
+    model::PrimitiveWet,
+)
+    large_scale_condensation!(column,scheme,
+        model.clausis_clapeyron,model.geometry,model.constants,model.time_stepping)
+end
+
+"""
+$(TYPEDSIGNATURES)
+Large-scale condensation for a `column` by relaxation back to a reference
+relative humidity if larger than that. Calculates the tendencies for
+specific humidity and temperature and integrates the large-scale
+precipitation vertically for output."""
+function large_scale_condensation!(
+    column::ColumnVariables{NF},
+    scheme::ImmediateCondensation,
+    clausius_clapeyron::AbstractClausiusClapeyron,
+    geometry::Geometry,
+    constants::DynamicsConstants,
+    time_stepping::TimeStepper,
+) where NF
+
+    (;temp, humid, pres) = column           # prognostic vars: specific humidity, pressure
+    (;temp_tend, humid_tend) = column       # tendencies to write into
+    (;sat_humid) = column                   # intermediate variable, calculated in thermodynamics!
+
+    # precompute scaling constant for precipitation output
+    pₛ = pres[end]                          # surface pressure
+    (;gravity, water_density) = constants
+    (;Δt_sec) = time_stepping
+    (;σ_levels_thick) = geometry
+    pₛΔt_gρ = pₛ*Δt_sec/gravity/water_density 
+
+    max_heating = scheme.max_heating/Δt_sec
+
+    @inbounds for k in eachindex(column)
+        if humid[k] > sat_humid[k]
+
+            # tendency for immediate humid = sat_humid
+            humid_tend_k = (sat_humid[k] - humid[k])/Δt_sec
+
+            # implicit correction, Frierson et al. 2006 eq. (21)
+            dqsat_dT = grad_saturation_humidity(clausius_clapeyron,temp[k],pres[k])
+            humid_tend_k /= 1 + scheme.latent_heat_cₚ[]*dqsat_dT
+
+            # latent heat release with maximum heating limiter for stability
+            temp_tend[k] += min(max_heating, -scheme.latent_heat_cₚ[] * humid_tend_k)
+
+            # If there is large-scale condensation at a level higher (i.e. smaller k) than
+            # the cloud-top previously diagnosed due to convection, then increase the cloud-top
+            column.cloud_top = min(column.cloud_top, k)             # Page 7 (last sentence)
+
+            # 2. Precipitation due to large-scale condensation [kg/m²/s] /ρ for [m/s]
+            # += for vertical integral
+            ΔpₖΔt_gρ = σ_levels_thick[k]*pₛΔt_gρ                    # Formula 4 *Δt for [m] of rain during Δt
+            column.precip_large_scale -= ΔpₖΔt_gρ * humid_tend_k    # Formula 25, unit [m]
+
             # only accumulate into humid_tend now to allow humid_tend != 0 before this scheme is called
             humid_tend[k] += humid_tend_k
         end

src/physics/surface_fluxes.jl

@@ -159,22 +159,22 @@ end
 # function barrier
 function evaporation!(  column::ColumnVariables,
                         model::PrimitiveWet)
-    evaporation!(column,model.evaporation,model.thermodynamics)
+    evaporation!(column,model.evaporation,model.clausis_clapeyron)
 end
 
 function evaporation!(  column::ColumnVariables{NF},
                         evaporation::SurfaceEvaporation,
-                        thermodynamics::Thermodynamics) where NF
+                        clausius_clapeyron::AbstractClausiusClapeyron) where NF
 
-    (;skin_temperature_sea, skin_temperature_land, pres, humid) = column
+    (;skin_temperature_sea, skin_temperature_land, pres) = column
     moisture_exchange_land = convert(NF,evaporation.moisture_exchange_land)
     moisture_exchange_sea  = convert(NF,evaporation.moisture_exchange_sea)
     α = column.soil_moisture_availability
 
     # SATURATION HUMIDITY OVER LAND AND OCEAN
     surface_pressure = pres[end]
-    sat_humid_land = saturation_humidity(skin_temperature_land,surface_pressure,thermodynamics)
-    sat_humid_sea = saturation_humidity(skin_temperature_sea,surface_pressure,thermodynamics)
+    sat_humid_land = saturation_humidity(skin_temperature_land,surface_pressure,clausius_clapeyron)
+    sat_humid_sea = saturation_humidity(skin_temperature_sea,surface_pressure,clausius_clapeyron)
 
     ρ = column.surface_air_density
     V₀ = column.surface_wind_speed

src/physics/tendencies.jl

@@ -13,6 +13,7 @@ function parameterization_tendencies!(
     model::PrimitiveEquation,
 )
 
+    # TODO move into shortwave radiation code
     cos_zenith!(time,model)
 
     G = model.geometry