diff --git a/docs/src/examples_3D.md b/docs/src/examples_3D.md index 66e111338..3fc556c61 100644 --- a/docs/src/examples_3D.md +++ b/docs/src/examples_3D.md @@ -212,6 +212,72 @@ And the comparison looks like ![Aquaplanet, no deep convection](aquaplanet_nodeepconvection.png) ![Aquaplanet, no convection](aquaplanet_noconvection.png) +## Large-scale vs convective precipitation + +```@example precipitation +using SpeedyWeather + +# components +spectral_grid = SpectralGrid(trunc=31, nlev=8) +large_scale_condensation = ImplicitCondensation(spectral_grid) +convection = SimplifiedBettsMiller(spectral_grid) + +# create model, initialize, run +model = PrimitiveWetModel(; spectral_grid, large_scale_condensation, convection) +simulation = initialize!(model) +model.feedback.verbose = false # hide +run!(simulation, period=Day(10)) +nothing # hide +``` + +We run the default `PrimitiveWetModel` with `ImplicitCondensation` as large-scale condensation +(see [Implicit large-scale condensation](@ref)) and the `SimplifiedBettsMiller` +for convection (see [Simplified Betts-Miller](@ref BettsMiller)). These schemes +have some additional parameters, we leave them as default for now, but you could +do `ImplicitCondensation(spectral_grid, relative_humidity_threshold = 0.8)` to +let it rain at 80% instead of 100% relative humidity. We now want to analyse +the precipitation that comes from these parameterizations + +```@example precipitation +using CairoMakie + +(; precip_large_scale, precip_convection) = simulation.diagnostic_variables.surface +m2mm = 1000 # convert from [m] to [mm] +heatmap(m2mm*precip_large_scale, title="Large-scale precipiation [mm]: Accumulated over 10 days", colormap=:dense) +save("large-scale_precipitation_acc.png", ans) # hide +nothing # hide +``` +![Large-scale precipitation](large-scale_precipitation_acc.png) + +Precipitation (both large-scale and convective) are written into the +`simulation.diagnostic_variables.surface` which, however, accumulate all precipitation +until netCDF output is written. As we didn't specify `output=true` here, they accumulated +precipitation over the last 10 days of the simulation which includes the initial +adjustment of humidity in the tropics (the large band of precipitation here). +So let us reset these accumulators and integrate for another 6 hours to get the +precipitation only in the period. + +```@example precipitation +# reset accumulators and simulate 6 hours +simulation.diagnostic_variables.surface.precip_large_scale .= 0 +simulation.diagnostic_variables.surface.precip_convection .= 0 +run!(simulation, period=Hour(6)) + +# visualise, precip_* arrays are flat copies, no need to read them out again! +m2mm_hr = (1000*Hour(1)/Hour(6)) # convert from [m] to [mm/hr] +heatmap(m2mm_hr*precip_large_scale, title="Large-scale precipiation [mm/hr]", colormap=:dense) +save("large-scale_precipitation.png", ans) # hide +heatmap(m2mm_hr*precip_convection, title="Convective precipiation [mm/hr]", colormap=:dense) +save("convective_precipitation.png", ans) # hide +nothing # hide +``` +![Large-scale precipitation](large-scale_precipitation.png) +![Convective precipitation](convective_precipitation.png) + +As the precipitation fields are accumulated meters over the integration period +(in the case of no output) we divide by 6 hours to get a precipitation rate ``[m/s]`` +but then multiply with 1 hour and 1000 to get the typical precipitation unit of ``[mm/hr]``. + ## References [^JW06]: Jablonowski, C. and Williamson, D.L. (2006), A baroclinic instability test case for atmospheric model dynamical cores. Q.J.R. Meteorol. Soc., 132: 2943-2975. DOI:[10.1256/qj.06.12](https://doi.org/10.1256/qj.06.12) diff --git a/src/dynamics/diagnostic_variables.jl b/src/dynamics/diagnostic_variables.jl index c35a4901e..4f4e7d162 100644 --- a/src/dynamics/diagnostic_variables.jl +++ b/src/dynamics/diagnostic_variables.jl @@ -35,7 +35,8 @@ Base.@kwdef struct GridVariables{NF<:AbstractFloat, Grid<:AbstractGrid{NF}} temp_grid ::Grid = zeros(Grid, nlat_half) # absolute temperature [K] temp_grid_prev ::Grid = zeros(Grid, nlat_half) # absolute temperature of previous time step [K] temp_virt_grid ::Grid = zeros(Grid, nlat_half) # virtual tempereature [K] - humid_grid ::Grid = zeros(Grid, nlat_half) # specific_humidity + humid_grid ::Grid = zeros(Grid, nlat_half) # specific_humidity [kg/kg] + humid_grid_prev ::Grid = zeros(Grid, nlat_half) # specific_humidity at previous time step u_grid ::Grid = zeros(Grid, nlat_half) # zonal velocity *coslat [m/s] v_grid ::Grid = zeros(Grid, nlat_half) # meridional velocity *coslat [m/s] u_grid_prev ::Grid = zeros(Grid, nlat_half) # zonal velocity *coslat of previous time step [m/s] diff --git a/src/dynamics/tendencies.jl b/src/dynamics/tendencies.jl index ad621a4c1..84e31c4a0 100644 --- a/src/dynamics/tendencies.jl +++ b/src/dynamics/tendencies.jl @@ -766,12 +766,14 @@ function SpeedyTransforms.gridded!( diagn::DiagnosticVariables, progn::PrognosticVariables, lf::Int, - model::ModelSetup) + model::ModelSetup; + kwargs... +) # all variables on layers for (progn_layer, diagn_layer) in zip(progn.layers, diagn.layers) progn_layer_lf = progn_layer.timesteps[lf] - gridded!(diagn_layer, progn_layer_lf, model) + gridded!(diagn_layer, progn_layer_lf, model; kwargs...) end # surface only for ShallowWaterModel or PrimitiveEquation @@ -788,9 +790,12 @@ $(TYPEDSIGNATURES) Propagate the spectral state of the prognostic variables `progn` to the diagnostic variables in `diagn` for the barotropic vorticity model. Updates grid vorticity, spectral stream function and spectral and grid velocities u, v.""" -function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer, - progn::PrognosticVariablesLayer, - model::Barotropic) +function SpeedyTransforms.gridded!( + diagn::DiagnosticVariablesLayer, + progn::PrognosticVariablesLayer, + model::Barotropic; + kwargs... +) (; vor_grid, u_grid, v_grid ) = diagn.grid_variables (; vor ) = progn # relative vorticity @@ -818,10 +823,12 @@ Propagate the spectral state of the prognostic variables `progn` to the diagnostic variables in `diagn` for the shallow water model. Updates grid vorticity, grid divergence, grid interface displacement (`pres_grid`) and the velocities u, v.""" -function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer, - progn::PrognosticVariablesLayer, - model::ShallowWater, # everything that's constant - ) +function SpeedyTransforms.gridded!( + diagn::DiagnosticVariablesLayer, + progn::PrognosticVariablesLayer, + model::ShallowWater; + kwargs... +) (; vor_grid, div_grid, u_grid, v_grid ) = diagn.grid_variables (; vor, div) = progn @@ -850,21 +857,27 @@ Propagate the spectral state of the prognostic variables `progn` to the diagnostic variables in `diagn` for primitive equation models. Updates grid vorticity, grid divergence, grid temperature, pressure (`pres_grid`) and the velocities u, v.""" -function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer, - progn::PrognosticVariablesLayer, - model::PrimitiveEquation) # everything that's constant +function SpeedyTransforms.gridded!( + diagn::DiagnosticVariablesLayer, + progn::PrognosticVariablesLayer, + model::PrimitiveEquation; + initialize=false +) (; vor_grid, div_grid, u_grid, v_grid ) = diagn.grid_variables (; temp_grid, humid_grid ) = diagn.grid_variables - (; temp_grid_prev, u_grid_prev, v_grid_prev) = diagn.grid_variables + (; temp_grid_prev, humid_grid_prev, u_grid_prev, v_grid_prev) = diagn.grid_variables (; vor, div, temp, humid) = progn U = diagn.dynamics_variables.a # reuse work arrays for velocities spectral V = diagn.dynamics_variables.b # U = u*coslat, V=v*coslat - # retain previous time step for vertical advection - @. temp_grid_prev = temp_grid - @. u_grid_prev = u_grid - @. v_grid_prev = v_grid + # retain previous time step for vertical advection and some parameterizations + if initialize == false # only store prev after initial step + @. temp_grid_prev = temp_grid # this is temperature anomaly wrt to implicit reference profile! + @. humid_grid_prev = humid_grid + @. u_grid_prev = u_grid + @. v_grid_prev = v_grid + end S = model.spectral_transform wet_core = model isa PrimitiveWet @@ -874,23 +887,31 @@ function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer, # V = v*coslat = coslat*∂ϕ/∂lat + ∂Ψ/dlon UV_from_vordiv!(U, V, vor, div, S) - gridded!(vor_grid, vor, S) # get vorticity on grid from spectral vor - gridded!(div_grid, div, S) # get divergence on grid from spectral div - gridded!(temp_grid, temp, S) # (absolute) temperature + gridded!(vor_grid, vor, S) # get vorticity on grid from spectral vor + gridded!(div_grid, div, S) # get divergence on grid from spectral div + gridded!(temp_grid, temp, S) # (absolute) temperature - if wet_core # specific humidity (wet core only) + if wet_core # specific humidity (wet core only) gridded!(humid_grid, humid, S) hole_filling!(humid_grid, model.hole_filling, model) # remove negative humidity end # include humidity effect into temp for everything stability-related temperature_average!(diagn, temp, S) - virtual_temperature!(diagn, temp, model) # temp = virt temp for dry core + virtual_temperature!(diagn, temp, model) # temp = virt temp for dry core # transform from U, V in spectral to u, v on grid (U, V = u, v*coslat) gridded!(u_grid, U, S, unscale_coslat=true) gridded!(v_grid, V, S, unscale_coslat=true) + if initialize # at initial step store prev <- current + # subtract the reference temperature profile as temp_grid is too after every time step + @. temp_grid_prev = temp_grid - model.implicit.temp_profile[diagn.k] + @. humid_grid_prev = humid_grid + @. u_grid_prev = u_grid + @. v_grid_prev = v_grid + end + return nothing end @@ -903,7 +924,6 @@ function temperature_average!( temp::LowerTriangularMatrix, S::SpectralTransform, ) - # average from l=m=0 harmonic divided by norm of the sphere diagn.temp_average[] = real(temp[1])/S.norm_sphere end \ No newline at end of file diff --git a/src/dynamics/time_integration.jl b/src/dynamics/time_integration.jl index 921638211..d1393a0b2 100644 --- a/src/dynamics/time_integration.jl +++ b/src/dynamics/time_integration.jl @@ -392,7 +392,7 @@ function time_stepping!( # propagate spectral state to grid variables for initial condition output (; output, feedback) = model lf = 1 # use first leapfrog index - gridded!(diagn, progn, lf, model) + gridded!(diagn, progn, lf, model, initialize=true) initialize!(progn.particles, progn, diagn, model.particle_advection) initialize!(output, feedback, time_stepping, clock, diagn, model) initialize!(model.callbacks, progn, diagn, model) diff --git a/src/physics/column_variables.jl b/src/physics/column_variables.jl index 684fd722f..28f4fff56 100644 --- a/src/physics/column_variables.jl +++ b/src/physics/column_variables.jl @@ -12,7 +12,8 @@ function get_column!( model.planet, model.orography, model.land_sea_mask, - model.albedo) + model.albedo, + model.implicit) end """ @@ -30,9 +31,11 @@ function get_column!( orography::AbstractOrography, land_sea_mask::AbstractLandSeaMask, albedo::AbstractAlbedo, + implicit::AbstractImplicit, ) (; σ_levels_full, ln_σ_levels_full) = geometry + (; temp_profile) = implicit # reference temperature on this layer @boundscheck C.nlev == D.nlev || throw(BoundsError) @@ -57,6 +60,12 @@ function get_column!( 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] end # TODO skin = surface approximation for now diff --git a/src/physics/define_column.jl b/src/physics/define_column.jl index a0c363a4b..92245b823 100644 --- a/src/physics/define_column.jl +++ b/src/physics/define_column.jl @@ -21,10 +21,16 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV orography::NF = 0 # orography height [m] # PROGNOSTIC VARIABLES - const u::Vector{NF} = zeros(NF, nlev) # zonal velocity [m/s] - const v::Vector{NF} = zeros(NF, nlev) # meridional velocity [m/s] - const temp::Vector{NF} = zeros(NF, nlev) # absolute temperature [K] - const humid::Vector{NF} = zeros(NF, nlev) # specific humidity [kg/kg] + const u::Vector{NF} = zeros(NF, nlev) # zonal velocity [m/s] + const v::Vector{NF} = zeros(NF, nlev) # meridional velocity [m/s] + 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)] diff --git a/src/physics/large_scale_condensation.jl b/src/physics/large_scale_condensation.jl index dc76f3da3..46480764b 100644 --- a/src/physics/large_scale_condensation.jl +++ b/src/physics/large_scale_condensation.jl @@ -11,11 +11,14 @@ export ImplicitCondensation Large scale condensation as with implicit precipitation. $(TYPEDFIELDS)""" Base.@kwdef struct ImplicitCondensation{NF<:AbstractFloat} <: AbstractCondensation - "Flux limiter for latent heat release [K] per timestep" - max_heating::NF = 0.2 + "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 = 1 - "Time scale in multiples of time step Δt" - time_scale::NF = 9 + "Time scale in multiples of time step Δt, the larger the less immediate" + time_scale::NF = 3 end ImplicitCondensation(SG::SpectralGrid; kwargs...) = ImplicitCondensation{SG.NF}(; kwargs...) @@ -58,18 +61,20 @@ relative humidity. Calculates the tendencies for specific humidity and temperature from latent heat release and integrates the large-scale precipitation vertically for output.""" function large_scale_condensation!( - column::ColumnVariables{NF}, + column::ColumnVariables, scheme::ImplicitCondensation, clausius_clapeyron::AbstractClausiusClapeyron, geometry::Geometry, planet::AbstractPlanet, atmosphere::AbstractAtmosphere, time_stepping::AbstractTimeStepper, -) 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! + (; 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 + (; 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 @@ -80,10 +85,10 @@ 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 = scheme.time_scale + (; time_scale, relative_humidity_threshold) = scheme @inbounds for k in eachindex(column) - if humid[k] > sat_humid[k] + if humid[k] > sat_humid[k]*relative_humidity_threshold # tendency for Implicit humid = sat_humid, divide by leapfrog time step below δq = sat_humid[k] - humid[k] diff --git a/src/physics/tendencies.jl b/src/physics/tendencies.jl index cb2267e5a..308819e61 100644 --- a/src/physics/tendencies.jl +++ b/src/physics/tendencies.jl @@ -36,6 +36,9 @@ function parameterization_tendencies!( end end +"""$(TYPEDSIGNATURES) +Calls for `column` one physics parameterization after another +and convert fluxes to tendencies.""" function parameterization_tendencies!( column::ColumnVariables, model::PrimitiveEquation diff --git a/src/physics/thermodynamics.jl b/src/physics/thermodynamics.jl index 500b06130..8e4c7da74 100644 --- a/src/physics/thermodynamics.jl +++ b/src/physics/thermodynamics.jl @@ -148,7 +148,11 @@ function saturation_humidity!( column::ColumnVariables, clausius_clapeyron::AbstractClausiusClapeyron, ) - (; sat_humid, pres, temp) = column + (; 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 for k in eachlayer(column) sat_humid[k] = saturation_humidity(temp[k], pres[k], clausius_clapeyron)