diff --git a/src/dynamics/atmosphere.jl b/src/dynamics/atmosphere.jl index f487a53fb..faccaa3e7 100644 --- a/src/dynamics/atmosphere.jl +++ b/src/dynamics/atmosphere.jl @@ -53,7 +53,10 @@ Base.@kwdef mutable struct EarthAtmosphere{NF<:AbstractFloat} <: AbstractAtmosph temp_ref::NF = 288 "reference moist-adiabatic temperature lapse rate [K/m]" - lapse_rate::NF = 5/1000 + moist_lapse_rate::NF = 5/1000 + + "reference dry-adiabatic temperature lapse rate [K/m]" + dry_lapse_rate::NF = 9.8/1000 "layer thickness for the shallow water model [m]" layer_thickness::NF = 8500 diff --git a/src/dynamics/initial_conditions.jl b/src/dynamics/initial_conditions.jl index 4e30edde1..19dd03e13 100644 --- a/src/dynamics/initial_conditions.jl +++ b/src/dynamics/initial_conditions.jl @@ -1,33 +1,55 @@ -abstract type InitialConditions end - -# EXPORT INITIAL CONDITIONS -export StartFromFile, - StartFromRest, - ZonalJet, - ZonalWind, - StartWithRandomVorticity +abstract type AbstractInitialConditions <: AbstractModelComponent end + +export InitialConditions +Base.@kwdef struct InitialConditions{V,P,T,H} <: AbstractInitialConditions + vordiv::V = ZeroInitially() + pres::P = ZeroInitially() + temp::T = ZeroInitially() + humid::H = ZeroInitially() +end -Base.@kwdef struct StartFromRest <: InitialConditions - pressure_on_orography::Bool = false +function initialize!( + progn::PrognosticVariables, + IC::InitialConditions, + model::ModelSetup +) + has(model,:vor) && initialize!(progn, IC.vordiv, model) + has(model,:pres) && initialize!(progn, IC.pres, model) + has(model,:temp) && initialize!(progn, IC.temp, model) + has(model,:humid) && initialize!(progn, IC.humid, model) end -function initialize!( progn::PrognosticVariables, - initial_conditions::StartFromRest, - model::ModelSetup) - return nothing # everything remains zero +InitialConditions(::Type{<:Barotropic}) = InitialConditions(;vordiv = StartWithRandomVorticity()) +InitialConditions(::Type{<:ShallowWater}) = InitialConditions(;vordiv = ZonalJet()) +function InitialConditions(::Type{<:PrimitiveDry}) + vordiv = ZonalWind() + pres = PressureOnOrography() + temp = JablonowskiTemperature() + return InitialConditions(;vordiv,pres,temp) end -function initialize!( progn::PrognosticVariables, - initial_conditions::StartFromRest, - model::PrimitiveEquation) - homogeneous_temperature!(progn, model) - initial_conditions.pressure_on_orography && pressure_on_orography!(progn, model) - # TODO initialise humidity +function InitialConditions(::Type{<:PrimitiveWet}) + vordiv = ZonalWind() + pres = PressureOnOrography() + temp = JablonowskiTemperature() + humid = ConstantRelativeHumidity() + return InitialConditions(;vordiv,pres,temp,humid) end +export ZeroInitially +struct ZeroInitially <: AbstractInitialConditions end +initialize!(::PrognosticVariables,::ZeroInitially,::ModelSetup) = nothing + +# to avoid a breaking change, like ZeroInitially +export StartFromRest +struct StartFromRest <: AbstractInitialConditions end +initialize!(::PrognosticVariables,::StartFromRest,::ModelSetup) = nothing + +export StartWithRandomVorticity + """Start with random vorticity as initial conditions $(TYPEDFIELDS)""" -Base.@kwdef mutable struct StartWithRandomVorticity <: InitialConditions +Base.@kwdef mutable struct StartWithRandomVorticity <: AbstractInitialConditions "Power of the spectral distribution k^power" power::Float64 = -3 @@ -57,12 +79,13 @@ function initialize!( progn::PrognosticVariables{NF}, end end +export ZonalJet + """ -$(TYPEDSIGNATURES) -Create a struct that contains all parameters for the Galewsky et al, 2004 zonal jet +A struct that contains all parameters for the Galewsky et al, 2004 zonal jet intitial conditions for the shallow water model. Default values as in Galewsky. $(TYPEDFIELDS)""" -Base.@kwdef mutable struct ZonalJet <: InitialConditions +Base.@kwdef mutable struct ZonalJet <: AbstractInitialConditions "jet latitude [˚N]" latitude::Float64 = 45 @@ -88,18 +111,12 @@ Base.@kwdef mutable struct ZonalJet <: InitialConditions perturb_height::Float64 = 120 end -function Base.show(io::IO, IC::InitialConditions) - println(io, "$(typeof(IC)) <: InitialConditions") - keys = propertynames(IC) - print_fields(io, IC, keys) -end - """ $(TYPEDSIGNATURES) Initial conditions from Galewsky, 2004, Tellus""" function initialize!( progn::PrognosticVariables, initial_conditions::ZonalJet, - model::ShallowWater) + model::ModelSetup) (; latitude, width, umax) = initial_conditions # for jet (; perturb_lat, perturb_lon, perturb_xwidth, # for perturbation @@ -175,17 +192,16 @@ function initialize!( progn::PrognosticVariables, spectral_truncation!(pres) end +export ZonalWind + """ $(TYPEDSIGNATURES) Create a struct that contains all parameters for the Jablonowski and Williamson, 2006 intitial conditions for the primitive equation model. Default values as in Jablonowski. $(TYPEDFIELDS)""" -Base.@kwdef struct ZonalWind <: InitialConditions +Base.@kwdef struct ZonalWind <: AbstractInitialConditions "conversion from σ to Jablonowski's ηᵥ-coordinates" η₀::Float64 = 0.252 - - "Sigma coordinates of the tropopause [1]" - σ_tropopause::Float64 = 0.2 "max amplitude of zonal wind [m/s]" u₀::Float64 = 35 @@ -202,20 +218,6 @@ Base.@kwdef struct ZonalWind <: InitialConditions "radius of Gaussian perturbation in units of Earth's radius [1]" perturb_radius::Float64 = 1/10 - - # TERMPERATURE - "reference dry-adiabtic lapse rate [K/m]" - lapse_rate::Float64 = 5/1000 - - "temperature difference used for stratospheric lapse rate [K], Jablonowski uses ΔT = 4.8e5 [K]" - ΔT::Float64 = 0 - - "minimum temperature [K] of profile" - Tmin::Float64 = 200 - - # PRESSURE - "initialize pressure given the `atmosphere.lapse_rate` on orography?" - pressure_on_orography::Bool = true end """ @@ -224,15 +226,11 @@ Initial conditions from Jablonowski and Williamson, 2006, QJR Meteorol. Soc""" function initialize!( progn::PrognosticVariables{NF}, initial_conditions::ZonalWind, model::PrimitiveEquation) where NF - - (; u₀, η₀, ΔT, Tmin, pressure_on_orography) = initial_conditions + + (; u₀, η₀) = initial_conditions (; perturb_lat, perturb_lon, perturb_uₚ, perturb_radius) = initial_conditions - (; lapse_rate, σ_tropopause) = initial_conditions - (; temp_ref, R_dry, pres_ref) = model.atmosphere - (; radius, Grid, nlat_half, nlev) = model.spectral_grid - (; rotation, gravity) = model.planet + (; radius, Grid, nlat_half) = model.spectral_grid (; σ_levels_full) = model.geometry - (; norm_sphere) = model.spectral_transform φ, λ = model.geometry.latds, model.geometry.londs S = model.spectral_transform @@ -283,11 +281,52 @@ function initialize!( progn::PrognosticVariables{NF}, spectral_truncation!(vor) spectral_truncation!(div) end +end - # TEMPERATURE - Tη = zero(σ_levels_full) +export JablonowskiTemperature + +""" +$(TYPEDSIGNATURES) +Create a struct that contains all parameters for the Jablonowski and Williamson, 2006 +intitial conditions for the primitive equation model. Default values as in Jablonowski. +$(TYPEDFIELDS)""" +Base.@kwdef struct JablonowskiTemperature <: AbstractInitialConditions + "conversion from σ to Jablonowski's ηᵥ-coordinates" + η₀::Float64 = 0.252 + + "Sigma coordinates of the tropopause [1]" + σ_tropopause::Float64 = 0.2 + + "max amplitude of zonal wind [m/s]" + u₀::Float64 = 35 + + "temperature difference used for stratospheric lapse rate [K], Jablonowski uses ΔT = 4.8e5 [K]" + ΔT::Float64 = 0 + + "minimum temperature [K] of profile" + Tmin::Float64 = 200 +end + +""" +$(TYPEDSIGNATURES) +Initial conditions from Jablonowski and Williamson, 2006, QJR Meteorol. Soc""" +function initialize!( progn::PrognosticVariables{NF}, + initial_conditions::JablonowskiTemperature, + model::ModelSetup) where NF + + (;u₀, η₀, ΔT, Tmin) = initial_conditions + (;σ_tropopause) = initial_conditions + lapse_rate = model.atmosphere.moist_lapse_rate + (;temp_ref, R_dry) = model.atmosphere + (;radius, Grid, nlat_half, nlev) = model.spectral_grid + (;rotation, gravity) = model.planet + (;σ_levels_full) = model.geometry + + φ, λ = model.geometry.latds, model.geometry.londs + S = model.spectral_transform # vertical profile + Tη = zero(σ_levels_full) for k in 1:nlev σ = σ_levels_full[k] Tη[k] = temp_ref*σ^(R_dry*lapse_rate/gravity) # Jablonowski and Williamson eq. 4 @@ -323,24 +362,16 @@ function initialize!( progn::PrognosticVariables{NF}, spectral!(temp, T, S) spectral_truncation!(temp) end - - # PRESSURE (constant everywhere) - lnp₀ = log(pres_ref) # logarithm of reference surface pressure [log(Pa)] - progn.surface.timesteps[1].pres[1] = norm_sphere*lnp₀ - - lnpₛ = ones(Grid{NF}, nlat_half) - lnpₛ .= pressure_on_orography ? pressure_on_orography!(progn, model) : lnp₀ - - # HUMIDITY - initialize_humidity!(progn, lnpₛ, model) end +export StartFromFile + """ Restart from a previous SpeedyWeather.jl simulation via the restart file restart.jld2 Applies interpolation in the horizontal but not in the vertical. restart.jld2 is identified by $(TYPEDFIELDS)""" -Base.@kwdef struct StartFromFile <: InitialConditions +Base.@kwdef struct StartFromFile <: AbstractInitialConditions "path for restart file" path::String = pwd() @@ -360,7 +391,11 @@ function initialize!( progn_new::PrognosticVariables, restart_file = jldopen(joinpath(path, string("run_", run_id_to_string(id)), "restart.jld2")) progn_old = restart_file["prognostic_variables"] - # version = restart_file["version"] # currently unused + version = restart_file["version"] + if version != pkgversion(SpeedyWeather) + @info "Restart file created with SpeedyWeather $version loaded"* + "but currently used is $(pkgversion(SpeedyWeather))" + end return copy!(progn_new, progn_old) end @@ -401,18 +436,27 @@ function homogeneous_temperature!( progn::PrognosticVariables, end end +export PressureOnOrography +struct PressureOnOrography <: AbstractInitialConditions end + """ $(TYPEDSIGNATURES) Initialize surface pressure on orography by integrating the hydrostatic equation with the reference temperature lapse rate.""" -function pressure_on_orography!(progn::PrognosticVariables, - model::PrimitiveEquation) - # temp_ref: Reference absolute T [K] at surface z = 0, constant lapse rate +function initialize!( progn::PrognosticVariables, + ::PressureOnOrography, + model::PrimitiveEquation) + + # temp_ref: Reference absolute T [K] at surface z = 0 # lapse_rate: Reference temperature lapse rate -dT/dz [K/m] # gravity: Gravitational acceleration [m/s^2] - # R: Specific gas constant for dry air [J/kg/K] + # R_dry: Specific gas constant for dry air [J/kg/K] # pres_ref: Reference surface pressure [hPa] - (; temp_ref, lapse_rate, pres_ref, R_dry ) = model.atmosphere + + (;atmosphere) = model + lapse_rate = model isa PrimitiveDry ? atmosphere.dry_lapse_rate : atmosphere.moist_lapse_rate + + (; temp_ref, pres_ref, R_dry ) = model.atmosphere (; gravity ) = model.planet (; orography ) = model.orography # orography on the grid @@ -429,52 +473,64 @@ function pressure_on_orography!(progn::PrognosticVariables, lnp = progn.surface.timesteps[1].pres spectral!(lnp, lnp_grid, model.spectral_transform) spectral_truncation!(lnp) # set lmax+1 row to zero - return lnp_grid # return grid for use in initialize_humidity! end -function initialize_humidity!( progn::PrognosticVariables, - pres_surf_grid::AbstractGrid, - model::PrimitiveDry) +export ConstantPressure +struct ConstantPressure <: AbstractInitialConditions end + +function initialize!( progn::PrognosticVariables, + ::ConstantPressure, + model::PrimitiveEquation) + (;pres_ref) = model.atmosphere + (;norm_sphere) = model.spectral_transform + + # logarithm of reference surface pressure [log(Pa)] + # set the l=m=0 mode, normalize correctly + progn.surface.timesteps[1].pres[1] = log(pres_ref) * norm_sphere + + # set other modes explicitly to zero + progn.surface.timesteps[1].pres[2:end] .= 0 return nothing end -function initialize_humidity!( - progn::PrognosticVariables{NF}, - pres_surf_grid::AbstractGrid, - model::PrimitiveWet -) where NF - - # TODO create a InitialHumidity <: InitialConditions to hold these parameters - relhumid_ref::NF = 0.7 - scale_height_humid::NF = 2.5 # scale height for specific humidity [km] - scale_height::NF = 7.5 # scale height for pressure [km] - scale_height_ratio = scale_height/scale_height_humid +export ConstantRelativeHumidity +Base.@kwdef struct ConstantRelativeHumidity <: AbstractInitialConditions + relhumid_ref::Float64 = 0.7 +end +function initialize!( + progn::PrognosticVariables, + IC::ConstantRelativeHumidity, + model::ModelSetup, +) + (; relhumid_ref) = IC (; nlev, σ_levels_full) = model.geometry + lnpₛ = progn.surface.timesteps[1].pres + pres_grid = gridded(lnpₛ, model.spectral_transform) + pres_grid .= exp.(pres_grid) - # 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) - for ij in eachgridpoint(humid_surf_grid) - q_sat = saturation_humidity(temp_grid[ij], exp(pres_surf_grid[ij]), model.clausius_clapeyron) - 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 levels above for k in 1:nlev - for lm in eachharmonic(humid_surf) - progn.layers[k].timesteps[1].humid[lm] = humid_surf[lm]*σ_levels_full[k]^scale_height_ratio + temp_grid = gridded(progn.layers[k].timesteps[1].temp, model.spectral_transform) + humid_grid = zero(temp_grid) + + for ij in eachgridpoint(humid_grid) + pₖ = σ_levels_full[k] * pres_grid[ij] + q_sat = saturation_humidity(temp_grid[ij], pₖ, model.clausius_clapeyron) + humid_grid[ij] = relhumid_ref*q_sat end + + (;humid) = progn.layers[k].timesteps[1] + spectral!(humid, humid_grid, model.spectral_transform) + spectral_truncation!(humid) end end +export RandomWaves + """Parameters for random initial conditions for the interface displacement η in the shallow water equations. $(TYPEDFIELDS)""" -Base.@kwdef struct RandomWaves <: InitialConditions +Base.@kwdef struct RandomWaves <: AbstractInitialConditions # random interface displacement field A::Float64 = 2000 # amplitude [m] lmin::Int64 = 10 # minimum wavenumber diff --git a/src/models/barotropic.jl b/src/models/barotropic.jl index 4767c4e61..5aeb21d9d 100644 --- a/src/models/barotropic.jl +++ b/src/models/barotropic.jl @@ -18,7 +18,7 @@ Base.@kwdef mutable struct BarotropicModel{ CO<:AbstractCoriolis, FR<:AbstractForcing, DR<:AbstractDrag, - IC<:InitialConditions, + IC<:AbstractInitialConditions, TS<:AbstractTimeStepper, ST<:SpectralTransform{NF}, IM<:AbstractImplicit, @@ -38,7 +38,7 @@ Base.@kwdef mutable struct BarotropicModel{ coriolis::CO = Coriolis(spectral_grid) forcing::FR = NoForcing() drag::DR = NoDrag() - initial_conditions::IC = StartWithRandomVorticity() + initial_conditions::IC = InitialConditions(Barotropic) # NUMERICS device_setup::DS = DeviceSetup(CPUDevice()) diff --git a/src/models/primitive_dry.jl b/src/models/primitive_dry.jl index 5e6ae8938..8d4efc618 100644 --- a/src/models/primitive_dry.jl +++ b/src/models/primitive_dry.jl @@ -14,7 +14,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{ GO<:AbstractGeopotential, OR<:AbstractOrography, AC<:AbstractAdiabaticConversion, - IC<:InitialConditions, + IC<:AbstractInitialConditions, LS<:AbstractLandSeaMask, OC<:AbstractOcean, LA<:AbstractLand, @@ -48,7 +48,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{ coriolis::CO = Coriolis(spectral_grid) geopotential::GO = Geopotential(spectral_grid) adiabatic_conversion::AC = AdiabaticConversion(spectral_grid) - initial_conditions::IC = ZonalWind() + initial_conditions::IC = InitialConditions(PrimitiveDry) # BOUNDARY CONDITIONS orography::OR = EarthOrography(spectral_grid) diff --git a/src/models/primitive_wet.jl b/src/models/primitive_wet.jl index 3d13304c5..bb4734e5f 100644 --- a/src/models/primitive_wet.jl +++ b/src/models/primitive_wet.jl @@ -14,7 +14,7 @@ Base.@kwdef mutable struct PrimitiveWetModel{ GO<:AbstractGeopotential, OR<:AbstractOrography, AC<:AbstractAdiabaticConversion, - IC<:InitialConditions, + IC<:AbstractInitialConditions, LS<:AbstractLandSeaMask, OC<:AbstractOcean, LA<:AbstractLand, @@ -55,7 +55,7 @@ Base.@kwdef mutable struct PrimitiveWetModel{ coriolis::CO = Coriolis(spectral_grid) geopotential::GO = Geopotential(spectral_grid) adiabatic_conversion::AC = AdiabaticConversion(spectral_grid) - initial_conditions::IC = ZonalWind() + initial_conditions::IC = InitialConditions(PrimitiveWet) # BOUNDARY CONDITIONS orography::OR = EarthOrography(spectral_grid) diff --git a/src/models/shallow_water.jl b/src/models/shallow_water.jl index e97503916..8f8ed9e9f 100644 --- a/src/models/shallow_water.jl +++ b/src/models/shallow_water.jl @@ -19,7 +19,7 @@ Base.@kwdef mutable struct ShallowWaterModel{ OR<:AbstractOrography, FR<:AbstractForcing, DR<:AbstractDrag, - IC<:InitialConditions, + IC<:AbstractInitialConditions, TS<:AbstractTimeStepper, ST<:SpectralTransform{NF}, IM<:AbstractImplicit, @@ -39,7 +39,7 @@ Base.@kwdef mutable struct ShallowWaterModel{ orography::OR = EarthOrography(spectral_grid) forcing::FR = NoForcing() drag::DR = NoDrag() - initial_conditions::IC = ZonalJet() + initial_conditions::IC = InitialConditions(ShallowWater) # NUMERICS time_stepping::TS = Leapfrog(spectral_grid) @@ -49,7 +49,7 @@ Base.@kwdef mutable struct ShallowWaterModel{ geometry::GE = Geometry(spectral_grid) # OUTPUT - output::OW = OutputWriter(spectral_grid, Barotropic) + output::OW = OutputWriter(spectral_grid, ShallowWater) callbacks::Dict{Symbol, AbstractCallback} = Dict{Symbol, AbstractCallback}() feedback::FB = Feedback() end