From 1a07c2a3f971e1bf53694bcd027dc1a06fc44c46 Mon Sep 17 00:00:00 2001 From: Milan Date: Mon, 23 Oct 2023 17:47:36 -0400 Subject: [PATCH] Overhauled initial conditions for humidity --- src/SpeedyWeather.jl | 2 +- src/dynamics/initial_conditions.jl | 19 ++++---- src/dynamics/models.jl | 4 +- src/physics/constants.jl | 61 ------------------------- src/physics/large_scale_condensation.jl | 11 +++-- src/physics/thermodynamics.jl | 6 +-- 6 files changed, 18 insertions(+), 85 deletions(-) delete mode 100644 src/physics/constants.jl diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl index 0cf01a018..0692ae031 100644 --- a/src/SpeedyWeather.jl +++ b/src/SpeedyWeather.jl @@ -96,7 +96,7 @@ export NoTemperatureRelaxation, JablonowskiRelaxation # EXPORT BOUNDARY LAYER SCHEMES -export NoBoundaryLayer, +export NoBoundaryLayerDrag, LinearDrag, QuadraticDrag diff --git a/src/dynamics/initial_conditions.jl b/src/dynamics/initial_conditions.jl index c05a03e9b..7a5bf423b 100644 --- a/src/dynamics/initial_conditions.jl +++ b/src/dynamics/initial_conditions.jl @@ -441,25 +441,22 @@ 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) @@ -467,7 +464,7 @@ function initialize_humidity!( progn::PrognosticVariables, # 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 diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl index 9f51332d5..1f2c7126a 100644 --- a/src/dynamics/models.jl +++ b/src/dynamics/models.jl @@ -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) @@ -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) diff --git a/src/physics/constants.jl b/src/physics/constants.jl deleted file mode 100644 index 5eff5857e..000000000 --- a/src/physics/constants.jl +++ /dev/null @@ -1,61 +0,0 @@ -struct ParameterizationConstants{NF<:AbstractFloat} <: AbstractParameterizationConstants{NF} - # BOUNDARY LAYER - drag_coefs::Vector{NF} # (inverse) drag time scale per layer - - # TEMPERATURE RELAXATION (Held and Suarez, 1996) - temp_relax_freq::Matrix{NF} # (inverse) relaxation time scale per layer and latitude - temp_equil_a::Vector{NF} # two terms to calculate equilibrium temperature as function - temp_equil_b::Vector{NF} # of latitude and pressure - - # TEMPERATURE RELAXATION (Jablonowski and Williamson, 2006) - temp_equil::Matrix{NF} - - # VERTICAL DIFFUSION - vert_diff_∇²_above::Vector{NF} - vert_diff_∇²_below::Vector{NF} #  - vert_diff_Δσ::Vector{NF} - - # RADIATION - fband::Matrix{NF} -end - -function Base.zeros(::Type{ParameterizationConstants}, - P::Parameters, - G::AbstractGeometry{NF}) where NF - (;nlev,nlat) = G - (;nband) = P - - drag_coefs = zeros(NF,nlev) # coefficient for boundary layer drag - - temp_relax_freq = zeros(NF,nlev,nlat) # (inverse) relaxation time scale per layer and latitude - temp_equil_a = zeros(NF,nlat) # terms to calculate equilibrium temperature as function - temp_equil_b = zeros(NF,nlat) # of latitude and pressure - - temp_equil = zeros(NF,nlev,nlat) - - vert_diff_∇²_above = zeros(NF,nlev-1) # defined on half levels, top, and bottom are =0 though - vert_diff_∇²_below = zeros(NF,nlev-1) # defined on half levels, top, and bottom are =0 though - vert_diff_Δσ = zeros(NF,nlev-1) # vertical gradient operator wrt σ coordinates - - fband = zeros(NF,400,nband) - - return ParameterizationConstants{NF}( drag_coefs, - temp_relax_freq, - temp_equil_a, - temp_equil_b, - temp_equil, - vert_diff_∇²_above, - vert_diff_∇²_below, - vert_diff_Δσ, - fband, - ) -end - -function ParameterizationConstants(P::Parameters,G::AbstractGeometry) - K = zeros(ParameterizationConstants,P,G) - initialize_boundary_layer!(K,P.boundary_layer,P,G) - initialize_temperature_relaxation!(K,P.temperature_relaxation,P,G) - initialize_vertical_diffusion!(K,P.vertical_diffusion,P,G) - initialize_longwave_radiation!(K,P) - return K -end diff --git a/src/physics/large_scale_condensation.jl b/src/physics/large_scale_condensation.jl index 4f4fee1cd..7af91eb11 100644 --- a/src/physics/large_scale_condensation.jl +++ b/src/physics/large_scale_condensation.jl @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 \ No newline at end of file diff --git a/src/physics/thermodynamics.jl b/src/physics/thermodynamics.jl index 771423a3c..0e0e6996b 100644 --- a/src/physics/thermodynamics.jl +++ b/src/physics/thermodynamics.jl @@ -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))