Skip to content

Commit

Permalink
Merge pull request #404 from SpeedyWeather/mk/humid
Browse files Browse the repository at this point in the history
Overhauled initial conditions for humidity
  • Loading branch information
milankl authored Oct 24, 2023
2 parents 59b674e + 1a07c2a commit 09d4a24
Show file tree
Hide file tree
Showing 6 changed files with 18 additions and 85 deletions.
2 changes: 1 addition & 1 deletion src/SpeedyWeather.jl
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@ export NoTemperatureRelaxation,
JablonowskiRelaxation

# EXPORT BOUNDARY LAYER SCHEMES
export NoBoundaryLayer,
export NoBoundaryLayerDrag,
LinearDrag,
QuadraticDrag

Expand Down
19 changes: 8 additions & 11 deletions src/dynamics/initial_conditions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -441,33 +441,30 @@ function initialize_humidity!( progn::PrognosticVariables,
pres_surf_grid::AbstractGrid,
model::PrimitiveWet)

# reference saturation water vapour pressure [Pa]
# relative humidity reference [1]
(;water_pres_ref, relhumid_ref, R_dry, R_vapour, pres_ref) = model.atmosphere
gas_ratio = R_dry/R_vapour
humid_ref = relhumid_ref*gas_ratio*water_pres_ref # reference specific humidity [Pa]
(;relhumid_ref) = model.atmosphere # relative humidity reference [1]

# ratio of scale heights [1], scale height [km], scale height for spec humidity [km]
(;scale_height, scale_height_humid, σ_tropopause) = model.atmosphere
scale_height_ratio = scale_height/scale_height_humid

(;nlev, σ_levels_full) = model.geometry
n_stratosphere_levels = findfirst->σ>=σ_tropopause,σ_levels_full)
n_stratosphere_levels = findlast->σ<=σ_tropopause,σ_levels_full)

# Specific humidity at the surface (grid space)
temp_grid = gridded(progn.layers[end].timesteps[1].temp,model.spectral_transform)
humid_surf_grid = zero(pres_surf_grid)
# @. humid_surf_grid = humid_ref*(exp(pres_surf_grid)/(pres_ref*100))^scale_height_ratio
q_ref = 1e-3 # kg/kg at the surface
@. humid_surf_grid .= q_ref
RingGrids.scale_coslat²!(humid_surf_grid)
for ij in eachgridpoint(humid_surf_grid)
q_sat = saturation_humidity(temp_grid[ij],exp(pres_surf_grid[ij]),model.thermodynamics)
humid_surf_grid[ij] = relhumid_ref*q_sat
end

humid_surf = spectral(humid_surf_grid,model.spectral_transform)
spectral_truncation!(humid_surf)

# Specific humidity at tropospheric levels (stratospheric humidity remains zero)
a = model.spectral_transform.norm_sphere
for k in n_stratosphere_levels+1:nlev
for lm in eachharmonic(humid_surf,progn.layers[1].timesteps[1].humid)
for lm in eachharmonic(humid_surf)
progn.layers[k].timesteps[1].humid[lm] = humid_surf[lm]*σ_levels_full[k]^scale_height_ratio
end
end
Expand Down
4 changes: 2 additions & 2 deletions src/dynamics/models.jl
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr

# PHYSICS/PARAMETERIZATIONS
physics::Bool = true
boundary_layer_drag::BoundaryLayerDrag{NF} = LinearDrag(spectral_grid)
boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid)
static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid)
surface_thermodynamics::AbstractSurfaceThermodynamics{NF} = SurfaceThermodynamicsConstant(spectral_grid)
Expand Down Expand Up @@ -248,7 +248,7 @@ Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr
# PHYSICS/PARAMETERIZATIONS
physics::Bool = true
thermodynamics::Thermodynamics{NF} = Thermodynamics(spectral_grid,atmosphere)
boundary_layer_drag::BoundaryLayerDrag{NF} = LinearDrag(spectral_grid)
boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid)
static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid)
large_scale_condensation::AbstractCondensation{NF} = SpeedyCondensation(spectral_grid)
Expand Down
61 changes: 0 additions & 61 deletions src/physics/constants.jl

This file was deleted.

11 changes: 6 additions & 5 deletions src/physics/large_scale_condensation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ function initialize!(scheme::SpeedyCondensation,model::PrimitiveEquation)
# the relative humidity threshold above which condensation occurs per layer
σₖ² = σ_levels_full[k]^2
relative_threshold[k] = threshold_max + threshold_range * (σₖ² - 1)
if k >= σ_boundary_layer
if σ_levels_full[k] >= σ_boundary_layer
relative_threshold[k] = max(relative_threshold[k], threshold_boundary_layer)
end

Expand Down Expand Up @@ -87,7 +87,7 @@ function large_scale_condensation!(
) where NF

(;relative_threshold,humid_tend_max) = scheme
time_scale = convert(NF,3600*scheme.time_scale)
time_scale = convert(NF,3600*scheme.time_scale) # [hrs] -> [s]
n_stratosphere_levels = scheme.n_stratosphere_levels[]

(;humid, pres) = column # prognostic variables: specific humidity, surface pressure
Expand All @@ -101,7 +101,7 @@ function large_scale_condensation!(
pₛΔt_gρ = pₛ*Δt_sec/gravity/water_density # precompute constant

(;σ_levels_thick) = geometry
latent_heat = convert(NF, atmosphere.latent_heat_condensation/atmosphere.cₚ)
latent_heat = convert(NF, atmosphere.latent_heat_condensation)

# 1. Tendencies of humidity and temperature due to large-scale condensation
@inbounds for k in n_stratosphere_levels+1:nlev # top to bottom, skip stratospheric levels
Expand All @@ -112,7 +112,8 @@ function large_scale_condensation!(
if humid[k] > humid_threshold
# accumulate in tendencies (nothing is added if humidity not above threshold)
humid_tend_k = -(humid[k] - humid_threshold) / time_scale # Formula 22
temp_tend[k] += -latent_heat * min(humid_tend_k, humid_tend_max[k]*pres[k]) # Formula 23
# temp_tend[k] += -latent_heat * min(humid_tend_k, humid_tend_max[k]*pres[k]) # Formula 23
temp_tend[k] += -latent_heat * humid_tend_k # Formula 23, without limiter

# 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
Expand All @@ -124,7 +125,7 @@ function large_scale_condensation!(
column.precip_large_scale += -ΔpₖΔt_gρ * humid_tend_k # Formula 25, unit [m]

# only write into humid_tend now to allow humid_tend != 0 before this scheme is called
humid_tend[k] += humid_tend_k
humid_tend[k] += humid_tend_k
end
end
end
6 changes: 1 addition & 5 deletions src/physics/thermodynamics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -88,20 +88,16 @@ function saturation_humidity!(
thermodynamics::Thermodynamics,
)
(;sat_humid, sat_vap_pres, pres, temp) = column
# (;e₀, T₀, C₁, C₂, T₁, T₂) = thermodynamics.magnus_coefs
(;mol_ratio) = thermodynamics # = mol_mass_vapour/mol_mass_dry_air = 0.622

for k in eachlayer(column)
# change coefficients for water (temp > T₀) or ice (else)
# C, T = temp[k] > T₀ ? (C₁, T₁) : (C₂, T₂)
# sat_vap_pres[k] = e₀ * exp(C * (temp[k] - T₀) / (temp[k] - T))
sat_vap_pres[k] = saturation_vapour_pressure(temp[k],thermodynamics.magnus_coefs)
# sat_humid[k] = mol_ratio*sat_vap_pres[k] / (pres[k] - (1-mol_ratio)*sat_vap_pres[k])
sat_humid[k] = saturation_humidity(sat_vap_pres[k],pres[k];mol_ratio)
end
end

function saturation_vapour_pressure(temp::NF,magnus_coefs::MagnusCoefs{NF}) where NF
# change coefficients for water (temp > T₀) or ice (else)
(;e₀, T₀, C₁, C₂, T₁, T₂) = magnus_coefs
C, T = temp > T₀ ? (C₁, T₁) : (C₂, T₂)
return e₀ * exp(C * (temp - T₀) / (temp - T))
Expand Down

0 comments on commit 09d4a24

Please sign in to comment.