From ec75a2e9206661bc669aa01fdf8428c9bc28e793 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= Date: Wed, 26 Jun 2024 12:46:15 +0100 Subject: [PATCH] All parameterizations at prev timestep, flux limiters removed, no shortwave (#555) * All parameterizations at prev timestep, flux limiters back to pre v0.10 * recalculate virtual temperature * Take all flux limiters out * Increase critical Richardson number to 10 * Remove all flux limiters * TransparentShortwave back in * surface heat flux max flux removed * NoShortwave as default for primitive models --- docs/src/surface_fluxes.md | 8 ------ src/dynamics/planet.jl | 6 ++--- src/dynamics/virtual_temperature.jl | 20 +++++++++++++++ src/models/primitive_dry.jl | 2 +- src/models/primitive_wet.jl | 2 +- src/physics/boundary_layer.jl | 2 +- src/physics/column_variables.jl | 22 +++++++---------- src/physics/define_column.jl | 6 ----- src/physics/large_scale_condensation.jl | 15 +++-------- src/physics/shortwave_radiation.jl | 5 +++- src/physics/surface_fluxes.jl | 33 ++++++------------------- src/physics/thermodynamics.jl | 9 +++---- src/physics/vertical_diffusion.jl | 6 ++--- 13 files changed, 58 insertions(+), 78 deletions(-) diff --git a/docs/src/surface_fluxes.md b/docs/src/surface_fluxes.md index 4f6eaa41d..0a3fbb7de 100644 --- a/docs/src/surface_fluxes.md +++ b/docs/src/surface_fluxes.md @@ -191,11 +191,6 @@ with ``\rho_s = \frac{p_s}{R_d T_N}`` the surface air density calculated from surface pressure ``p_s`` and lowermost layer temperature ``T_N``. Better would be to extrapolate ``T_N`` to ``T_s`` a surface air temperature assuming adiabatic descent but that is currently not implemented. -In practice we use a flux limiter for numerical stability which limits the magnitude -(preserving the sign). Choosing ``F_{uv, max}^\uparrow = 0.5 Pa`` would mean -that the drag for winds faster than about ``33~m/s`` (typical ``C = 5\cdot 10^{-4}``) -does not further increase. This can help to prevent oscillations drag terms can produce -for sudden strong wind gusts. ## Surface heat fluxes @@ -212,9 +207,6 @@ land surface layer or the mixed layer in the ocean. We then compute F_T^\uparrow = \rho_s C V_0 c_p (T_{skin} - T_s) ``` -and apply a similar flux limiter as for the momentum flux to prevent a sudden strong heating -or cooling. For ``F_{T, max}^\uparrow = 100~W/m^2`` - ## Surface evaporation The surface moisture flux, i.e. evaporation of soil moisture over land and evaporation of diff --git a/src/dynamics/planet.jl b/src/dynamics/planet.jl index 45e34e796..f11874f1e 100644 --- a/src/dynamics/planet.jl +++ b/src/dynamics/planet.jl @@ -1,6 +1,6 @@ abstract type AbstractPlanet <: AbstractModelComponent end -const DEFAULT_ROTATION = 7.29e-5 # default angular frequency of Earth's rotation [1/s] +const DEFAULT_ROTATION = 7.29e-5 # default angular frequency of Earth's rotation [1/s] const DEFAULT_GRAVITY = 9.81 # default gravitational acceleration on Earth [m/s²] export Earth @@ -12,7 +12,7 @@ characteristics. Note that `radius` is not part of it as this should be chosen in `SpectralGrid`. Keyword arguments are $(TYPEDFIELDS) """ -Base.@kwdef mutable struct Earth{NF<:AbstractFloat} <: AbstractPlanet +@kwdef mutable struct Earth{NF<:AbstractFloat} <: AbstractPlanet "angular frequency of Earth's rotation [rad/s]" rotation::NF = DEFAULT_ROTATION @@ -39,7 +39,7 @@ Base.@kwdef mutable struct Earth{NF<:AbstractFloat} <: AbstractPlanet axial_tilt::NF = 23.4 "Total solar irradiance at the distance of 1 AU [W/m²]" - solar_constant::NF = 1365/4 # for testing + solar_constant::NF = 1365 end Earth(SG::SpectralGrid; kwargs...) = Earth{SG.NF}(; kwargs...) diff --git a/src/dynamics/virtual_temperature.jl b/src/dynamics/virtual_temperature.jl index 5b7768893..a99c71410 100644 --- a/src/dynamics/virtual_temperature.jl +++ b/src/dynamics/virtual_temperature.jl @@ -50,6 +50,26 @@ function virtual_temperature!( diagn::DiagnosticVariablesLayer, copyto!(temp_virt_grid, temp_grid) end +function virtual_temperature!( + column::ColumnVariables, + model::PrimitiveEquation, +) + (; temp, temp_virt, humid) = column + μ = model.atmosphere.μ_virt_temp + + @. temp_virt = temp*(1 + μ*humid) + return nothing +end + +function virtual_temperature!( + column::ColumnVariables, + model::PrimitiveDry, +) + (; temp, temp_virt) = column + @. temp_virt = temp # temp = temp_virt for PrimitiveDry + return nothing +end + """ $(TYPEDSIGNATURES) Linear virtual temperature for `model::PrimitiveDry`: Just copy over diff --git a/src/models/primitive_dry.jl b/src/models/primitive_dry.jl index 4e87258f1..b9ebece98 100644 --- a/src/models/primitive_dry.jl +++ b/src/models/primitive_dry.jl @@ -76,7 +76,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{ surface_wind::SUW = SurfaceWind(spectral_grid) surface_heat_flux::SH = SurfaceHeatFlux(spectral_grid) convection::CV = DryBettsMiller(spectral_grid) - shortwave_radiation::SW = TransparentShortwave(spectral_grid) + shortwave_radiation::SW = NoShortwave(spectral_grid) longwave_radiation::LW = JeevanjeeRadiation(spectral_grid) # NUMERICS diff --git a/src/models/primitive_wet.jl b/src/models/primitive_wet.jl index 7ff800f63..010028419 100644 --- a/src/models/primitive_wet.jl +++ b/src/models/primitive_wet.jl @@ -87,7 +87,7 @@ Base.@kwdef mutable struct PrimitiveWetModel{ surface_evaporation::EV = SurfaceEvaporation(spectral_grid) large_scale_condensation::LSC = ImplicitCondensation(spectral_grid) convection::CV = SimplifiedBettsMiller(spectral_grid) - shortwave_radiation::SW = TransparentShortwave(spectral_grid) + shortwave_radiation::SW = NoShortwave(spectral_grid) longwave_radiation::LW = JeevanjeeRadiation(spectral_grid) # NUMERICS diff --git a/src/physics/boundary_layer.jl b/src/physics/boundary_layer.jl index 050bbc3b3..41d6202b5 100644 --- a/src/physics/boundary_layer.jl +++ b/src/physics/boundary_layer.jl @@ -102,7 +102,7 @@ Base.@kwdef struct BulkRichardsonDrag{NF} <: AbstractBoundaryLayer z₀::NF = 3.21e-5 "Critical Richardson number for stable mixing cutoff [1]" - Ri_c::NF = 1 + Ri_c::NF = 10 "Maximum drag coefficient, κ²/log(zₐ/z₀) for zₐ from reference temperature" drag_max::Base.RefValue{NF} = Ref(zero(NF)) diff --git a/src/physics/column_variables.jl b/src/physics/column_variables.jl index 28f4fff56..a79c3460e 100644 --- a/src/physics/column_variables.jl +++ b/src/physics/column_variables.jl @@ -35,7 +35,7 @@ function get_column!( ) (; σ_levels_full, ln_σ_levels_full) = geometry - (; temp_profile) = implicit # reference temperature on this layer + (; temp_profile) = implicit # reference temperature @boundscheck C.nlev == D.nlev || throw(BoundsError) @@ -54,18 +54,14 @@ function get_column!( C.pres[1:end-1] .= σ_levels_full.*pₛ # pressure on every level p = σ*pₛ C.pres[end] = pₛ # last element is surface pressure pₛ - @inbounds for (k, layer) = enumerate(D.layers) # read out prognostic variables on grid - C.u[k] = layer.grid_variables.u_grid[ij] - C.v[k] = layer.grid_variables.v_grid[ij] - C.temp[k] = layer.grid_variables.temp_grid[ij] - C.temp_virt[k] = layer.grid_variables.temp_virt_grid[ij] # actually diagnostic - C.humid[k] = layer.grid_variables.humid_grid[ij] - - # and at previous time step, add temp reference profile back in as temp_grid_prev is anomaly - C.u_prev[k] = layer.grid_variables.u_grid_prev[ij] - C.v_prev[k] = layer.grid_variables.v_grid_prev[ij] - C.temp_prev[k] = layer.grid_variables.temp_grid_prev[ij] + temp_profile[k] - C.humid_prev[k] = layer.grid_variables.humid_grid_prev[ij] + @inbounds for (k, layer) = enumerate(D.layers) + # read out prognostic variables on grid at previous time step + C.u[k] = layer.grid_variables.u_grid_prev[ij] + C.v[k] = layer.grid_variables.v_grid_prev[ij] + + # add temp reference profile back in as temp_grid_prev is anomaly + C.temp[k] = layer.grid_variables.temp_grid_prev[ij] + temp_profile[k] + C.humid[k] = layer.grid_variables.humid_grid_prev[ij] end # TODO skin = surface approximation for now diff --git a/src/physics/define_column.jl b/src/physics/define_column.jl index 92245b823..a93869041 100644 --- a/src/physics/define_column.jl +++ b/src/physics/define_column.jl @@ -26,12 +26,6 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV const temp::Vector{NF} = zeros(NF, nlev) # absolute temperature [K] const humid::Vector{NF} = zeros(NF, nlev) # specific humidity [kg/kg] - # PROGNOSTIC VARIABLES at previous time step - const u_prev::Vector{NF} = zeros(NF, nlev) # zonal velocity [m/s] - const v_prev::Vector{NF} = zeros(NF, nlev) # meridional velocity [m/s] - const temp_prev::Vector{NF} = zeros(NF, nlev) # absolute temperature [K] - const humid_prev::Vector{NF} = zeros(NF, nlev) # specific humidity [kg/kg] - # (log) pressure per layer, surface is prognostic, last element here, but precompute other layers too const ln_pres::Vector{NF} = zeros(NF, nlev+1) # logarithm of pressure [log(Pa)] const pres::Vector{NF} = zeros(NF, nlev+1) # pressure [Pa] diff --git a/src/physics/large_scale_condensation.jl b/src/physics/large_scale_condensation.jl index 73895b852..8829478c1 100644 --- a/src/physics/large_scale_condensation.jl +++ b/src/physics/large_scale_condensation.jl @@ -10,12 +10,9 @@ export ImplicitCondensation """ Large scale condensation as with implicit precipitation. $(TYPEDFIELDS)""" -Base.@kwdef struct ImplicitCondensation{NF<:AbstractFloat} <: AbstractCondensation +@kwdef struct ImplicitCondensation{NF<:AbstractFloat} <: AbstractCondensation "Relative humidity threshold [1 = 100%] to trigger condensation" relative_humidity_threshold::NF = 1 - - "Flux limiter for latent heat release [W/m²] per timestep" - max_heating::NF = 10 # currently not needed? maybe too high? "Time scale in multiples of time step Δt, the larger the less immediate" time_scale::NF = 3 @@ -70,9 +67,7 @@ function large_scale_condensation!( time_stepping::AbstractTimeStepper, ) - (; pres) = column # prognostic vars: pressure - temp = column.temp_prev # but use temp, humid from - humid = column.humid_prev # from previous time step for numerical stability + (; pres, temp, humid) = column # prognostic vars (from previous time step for numerical stability) (; temp_tend, humid_tend) = column # tendencies to write into (; sat_humid) = column # intermediate variable, calculated in thermodynamics! @@ -84,7 +79,6 @@ function large_scale_condensation!( (; Lᵥ, cₚ, Lᵥ_Rᵥ) = clausius_clapeyron Lᵥ_cₚ = Lᵥ/cₚ # latent heat of vaporization over heat capacity - max_heating = scheme.max_heating/Δt_sec (; time_scale, relative_humidity_threshold) = scheme @inbounds for k in eachindex(column) @@ -97,9 +91,8 @@ function large_scale_condensation!( dqsat_dT = sat_humid[k] * Lᵥ_Rᵥ/temp[k]^2 # derivative of qsat wrt to temp δq /= ((1 + Lᵥ_cₚ*dqsat_dT) * time_scale*Δt_sec) - # latent heat release with maximum heating limiter for stability - δT = min(max_heating, -Lᵥ_cₚ * δq) - δq = -δT/Lᵥ_cₚ # also limit drying for enthalpy conservation + # latent heat release for enthalpy conservation + δT = -Lᵥ_cₚ * δq # 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 diff --git a/src/physics/shortwave_radiation.jl b/src/physics/shortwave_radiation.jl index 6293d719c..3988cd424 100644 --- a/src/physics/shortwave_radiation.jl +++ b/src/physics/shortwave_radiation.jl @@ -34,5 +34,8 @@ function shortwave_radiation!( ) (; cos_zenith, albedo) = column (; solar_constant) = planet - column.flux_temp_upward[end] += (1-albedo) * solar_constant * cos_zenith + + # reduce strength by 1/4 as this currently just hits the air temperature in lowermost + # layer directly and not through a skin temperature + column.flux_temp_upward[end] += ((1-albedo) * solar_constant * cos_zenith)/4 end \ No newline at end of file diff --git a/src/physics/surface_fluxes.jl b/src/physics/surface_fluxes.jl index 634d1ef7d..5809013c6 100644 --- a/src/physics/surface_fluxes.jl +++ b/src/physics/surface_fluxes.jl @@ -29,8 +29,8 @@ function surface_thermodynamics!( column::ColumnVariables, # surface value is same as lowest model level, use previous time step # for numerical stability - column.surface_temp = column.temp_prev[end] # todo use constant POTENTIAL temperature - column.surface_humid = column.humid_prev[end] # humidity at surface is the same as + column.surface_temp = column.temp[end] # todo use constant POTENTIAL temperature + column.surface_humid = column.humid[end] # humidity at surface is the same as # surface air density via virtual temperature (; R_dry) = model.atmosphere @@ -44,7 +44,7 @@ function surface_thermodynamics!( column::ColumnVariables, (; R_dry) = model.atmosphere # surface value is same as lowest model level, but use previous # time step for numerical stability - column.surface_temp = column.temp_prev[end] # todo use constant POTENTIAL temperature + column.surface_temp = column.temp[end] # todo use constant POTENTIAL temperature column.surface_air_density = column.pres[end]/(R_dry*column.surface_temp) end @@ -71,9 +71,6 @@ Base.@kwdef struct SurfaceWind{NF<:AbstractFloat} <: AbstractSurfaceWind "Otherwise, Drag coefficient over sea [1]" drag_sea::NF = 1.8e-3 - - "Flux limiter to cap the max of surface momentum fluxes [kg/m/s²]" - max_flux::NF = 10 # currently not needed? maybe too high? end SurfaceWind(SG::SpectralGrid; kwargs...) = SurfaceWind{SG.NF}(; kwargs...) @@ -85,11 +82,10 @@ function surface_wind_stress!( column::ColumnVariables, (; land_fraction) = column (; f_wind, V_gust, drag_land, drag_sea) = surface_wind - (; max_flux) = surface_wind # SPEEDY documentation eq. 49, but use previous time step for numerical stability - column.surface_u = f_wind*column.u_prev[end] - column.surface_v = f_wind*column.v_prev[end] + column.surface_u = f_wind*column.u[end] + column.surface_v = f_wind*column.v[end] (; surface_u, surface_v) = column # SPEEDY documentation eq. 50 @@ -105,15 +101,9 @@ function surface_wind_stress!( column::ColumnVariables, V₀ = column.surface_wind_speed drag = land_fraction*drag_land + (1-land_fraction)*drag_sea - # add flux limiter to avoid heavy drag in (initial) shock - u_flux = ρ*drag*V₀*surface_u - v_flux = ρ*drag*V₀*surface_v - column.flux_u_upward[end] -= clamp(u_flux, -max_flux, max_flux) - column.flux_v_upward[end] -= clamp(v_flux, -max_flux, max_flux) - # SPEEDY documentation eq. 52, 53, accumulate fluxes with += - # column.flux_u_upward[end] -= ρ*drag*V₀*surface_u - # column.flux_v_upward[end] -= ρ*drag*V₀*surface_v + column.flux_u_upward[end] -= ρ*drag*V₀*surface_u + column.flux_v_upward[end] -= ρ*drag*V₀*surface_v return nothing end @@ -136,9 +126,6 @@ Base.@kwdef struct SurfaceHeatFlux{NF<:AbstractFloat} <: AbstractSurfaceHeatFlux "Otherwise, use the following drag coefficient for heat fluxes over sea" heat_exchange_sea::NF = 0.9e-3 - - "Flux limiter for surface heat fluxes [W/m²]" - max_flux::NF = 100 # currently not needed? maybe too high? end SurfaceHeatFlux(SG::SpectralGrid; kwargs...) = SurfaceHeatFlux{SG.NF}(; kwargs...) @@ -150,7 +137,7 @@ function surface_heat_flux!( model::PrimitiveEquation, ) cₚ = model.atmosphere.heat_capacity - (; heat_exchange_land, heat_exchange_sea, max_flux) = heat_flux + (; heat_exchange_land, heat_exchange_sea) = heat_flux ρ = column.surface_air_density V₀ = column.surface_wind_speed @@ -168,10 +155,6 @@ function surface_heat_flux!( flux_land = ρ*drag_land*V₀*cₚ*(T_skin_land - T) flux_sea = ρ*drag_sea*V₀*cₚ*(T_skin_sea - T) - # flux limiter - flux_land = clamp(flux_land, -max_flux, max_flux) - flux_sea = clamp(flux_sea, -max_flux, max_flux) - # mix fluxes for fractional land-sea mask land_available = isfinite(T_skin_land) sea_available = isfinite(T_skin_sea) diff --git a/src/physics/thermodynamics.jl b/src/physics/thermodynamics.jl index 8e4c7da74..4753694a8 100644 --- a/src/physics/thermodynamics.jl +++ b/src/physics/thermodynamics.jl @@ -121,6 +121,7 @@ Calculate geopotentiala and dry static energy for the primitive equation model." function get_thermodynamics!(column::ColumnVariables, model::PrimitiveEquation) geopotential!(column.geopot, column.temp, model.geopotential, column.surface_geopotential) dry_static_energy!(column, model.atmosphere) + virtual_temperature!(column, model) end """ @@ -148,13 +149,11 @@ function saturation_humidity!( column::ColumnVariables, clausius_clapeyron::AbstractClausiusClapeyron, ) - (; sat_humid, pres) = column - # use previous time step for temperature for stability of large-scale condensation # TODO also use previous pressure, but sat_humid is only weakly dependent on it, skip for now - temp = column.temp_prev + (; sat_humid, pres, temp) = column - for k in eachlayer(column) + @inbounds for k in eachlayer(column) sat_humid[k] = saturation_humidity(temp[k], pres[k], clausius_clapeyron) end end @@ -175,7 +174,7 @@ function moist_static_energy!( (; sat_moist_static_energy, moist_static_energy, dry_static_energy) = column (; humid, sat_humid) = column - for k in eachlayer(column) + @inbounds for k in eachlayer(column) moist_static_energy[k] = dry_static_energy[k] + Lᵥ * humid[k] sat_moist_static_energy[k] = dry_static_energy[k] + Lᵥ * sat_humid[k] end diff --git a/src/physics/vertical_diffusion.jl b/src/physics/vertical_diffusion.jl index f392af484..ecb2e7aee 100644 --- a/src/physics/vertical_diffusion.jl +++ b/src/physics/vertical_diffusion.jl @@ -16,7 +16,7 @@ initialize!(::NoVerticalDiffusion,::PrimitiveEquation) = nothing vertical_diffusion!(::ColumnVariables, ::NoVerticalDiffusion, ::PrimitiveEquation) = nothing export BulkRichardsonDiffusion -Base.@kwdef struct BulkRichardsonDiffusion{NF} <: AbstractVerticalDiffusion +@kwdef struct BulkRichardsonDiffusion{NF} <: AbstractVerticalDiffusion nlev::Int "[OPTION] von Kármán constant [1]" @@ -26,7 +26,7 @@ Base.@kwdef struct BulkRichardsonDiffusion{NF} <: AbstractVerticalDiffusion z₀::NF = 3.21e-5 "[OPTION] Critical Richardson number for stable mixing cutoff [1]" - Ri_c::NF = 1 + Ri_c::NF = 10 "[OPTION] Fraction of surface boundary layer" fb::NF = 0.1 @@ -216,7 +216,7 @@ function bulk_richardson!( ) cₚ = atmosphere.heat_capacity (; u, v, geopot, temp_virt, nlev) = column - bulk_richardson = column.a # reuse work array + bulk_richardson = column.d # reuse work array # surface layer V² = u[nlev]^2 + v[nlev]^2