From c8a38e0de1b5b996ac0f861245258b57002a209e Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 20:11:32 -0400 Subject: [PATCH 01/14] Split into 2D/3D examples --- docs/make.jl | 3 +- docs/src/convection.md | 2 +- docs/src/{examples.md => examples_2D.md} | 105 +-------------- docs/src/examples_3D.md | 160 +++++++++++++++++++++++ docs/src/forcing_drag.md | 2 +- docs/src/how_to_run_speedy.md | 4 +- docs/src/installation.md | 2 +- 7 files changed, 174 insertions(+), 104 deletions(-) rename docs/src/{examples.md => examples_2D.md} (75%) create mode 100644 docs/src/examples_3D.md diff --git a/docs/make.jl b/docs/make.jl index 36310d006..edec31989 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -27,7 +27,8 @@ makedocs( ], "Running SpeedyWeather" => [ "How to run SpeedyWeather"=>"how_to_run_speedy.md", - "Examples"=>"examples.md", + "Examples 2D"=>"examples_2D.md", + "Examples 3D"=>"examples_3D.md", "Analysis"=>"analysis.md", "Tree structure"=>"structure.md", "Particle advection"=>"particles.md", diff --git a/docs/src/convection.md b/docs/src/convection.md index c24f0d482..eb452fbf7 100644 --- a/docs/src/convection.md +++ b/docs/src/convection.md @@ -15,7 +15,7 @@ precipitating as the condensed humidity forms cloud droplets that eventually fall down as convective precipitation. See also [Large-scale condensation](@ref) in comparison. -## Simplified Betts-Miller convective adjustment +## [Simplified Betts-Miller convective adjustment](@ref BettsMiller) We follow the simplification of the Betts-Miller convection scheme [^Betts1986][^BettsMiller1986] as studied by Frierson, 2007 [^Frierson2007]. diff --git a/docs/src/examples.md b/docs/src/examples_2D.md similarity index 75% rename from docs/src/examples.md rename to docs/src/examples_2D.md index 4b90d3796..3b968dd0d 100644 --- a/docs/src/examples.md +++ b/docs/src/examples_2D.md @@ -1,8 +1,10 @@ -# Examples +# [Examples 2D](@id Examples) -The following is a collection of model setups, starting with an easy setup +The following is a collection of example model setups, starting with an easy setup of the [Barotropic vorticity equation](@ref barotropic_vorticity_model) and -continuing with more complicated setups. +continuing with the [shallow water model](@ref shallow_water_model). + +See also [Examples 3D](@ref) for examples with the primitive equation models. ## 2D turbulence on a non-rotating sphere @@ -304,7 +306,7 @@ Hb = model.orography.orography η = simulation.diagnostic_variables.surface.pres_grid h = @. η + H - Hb # @. to broadcast grid + scalar - grid -heatmap(h, title="Dynamic layer thickness h") +heatmap(h, title="Dynamic layer thickness h", colormap=:oslo) save("gravity_waves.png", ans) # hide nothing # hide ``` @@ -312,99 +314,6 @@ nothing # hide Can you spot the Himalayas or the Andes? -## Jablonowski-Williamson baroclinic wave - -```@example jablonowski -using SpeedyWeather -spectral_grid = SpectralGrid(trunc=31, nlev=8, Grid=FullGaussianGrid, dealiasing=3) - -orography = ZonalRidge(spectral_grid) -initial_conditions = InitialConditions( - vordiv = ZonalWind(), - temp = JablonowskiTemperature(), - pres = ZeroInitially()) - -model = PrimitiveDryModel(; spectral_grid, orography, initial_conditions, physics=false) -simulation = initialize!(model) -model.feedback.verbose = false # hide -run!(simulation, period=Day(9)) -nothing # hide -``` - -The Jablonowski-Williamson baroclinic wave test case[^JW06] using the -[Primitive equation model](@ref primitive_equation_model) particularly the dry model, -as we switch off all physics with `physics=false`. -We want to use 8 vertical levels, and a lower resolution of T31 on a -[full Gaussian grid](@ref FullGaussianGrid). -The Jablonowski-Williamson initial conditions are `ZonalWind` for vorticity and divergence -(curl and divergence of ``u, v``), `JablonowskiTemperature` for temperature and -`ZeroInitially` for pressure. The orography is just a `ZonalRidge`. -There is no forcing and the initial conditions are baroclinically unstable which kicks -off a wave propagating eastward. This wave becomes obvious when visualised with - -```@example jablonowski -using CairoMakie - -vor = simulation.diagnostic_variables.layers[end].grid_variables.vor_grid -heatmap(vor, title="Surface relative vorticity") -save("jablonowski.png", ans) # hide -nothing # hide -``` -![Jablonowski pyplot](jablonowski.png) - -## Aquaplanet - -```@example aquaplanet -using SpeedyWeather - -# components -spectral_grid = SpectralGrid(trunc=31, nlev=5) -ocean = AquaPlanet(spectral_grid, temp_equator=302, temp_poles=273) -land_sea_mask = AquaPlanetMask(spectral_grid) -orography = NoOrography(spectral_grid) - -# create model, initialize, run -model = PrimitiveWetModel(; spectral_grid, ocean, land_sea_mask, orography) -simulation = initialize!(model) -model.feedback.verbose = false # hide -run!(simulation, period=Day(50)) -nothing # hide -``` - -Here we have defined an aquaplanet simulation by -- creating an `ocean::AquaPlanet`. This will use constant sea surface temperatures that only vary with latitude. -- creating a `land_sea_mask::AquaPlanetMask` this will use a land-sea mask with `false`=ocean everywhere. -- creating an `orography::NoOrography` which will have no orography and zero surface geopotential. - -All passed on to the model constructor for a `PrimitiveWetModel`, we have now a model with humidity -and physics parameterization as they are defined by default (typing `model` will give you an overview -of its components). We could have change the `model.land` and `model.vegetation` components too, -but given the land-sea masks masks those contributions to the surface fluxes anyway, this is not -necessary. Note that neither sea surface temperature, land-sea mask -or orography have to agree. It is possible to have an ocean on top of a mountain. -For an ocean grid-cell that is (partially) masked by the land-sea mask, its value will -be (fractionally) ignored in the calculation of surface fluxes (potentially leading -to a zero flux depending on land surface temperatures). - -Now with the following we visualize the surface humidity after the 50 days of -simulation. We use 50 days as without mountains it takes longer for the initial conditions to -become unstable. The surface humidity shows small-scale patches in the tropics, which is a result -of the convection scheme, causing updrafts and downdrafts in both humidity and temperature. - -```@example aquaplanet -using CairoMakie - -humid = simulation.diagnostic_variables.layers[end].grid_variables.humid_grid -heatmap(humid, title="Surface specific humidity [kg/kg]") - -save("aquaplanet.png", ans) # hide -nothing # hide -``` -![Aquaplanet pyplot](aquaplanet.png) - - - ## References -[^G04]: Galewsky, Scott, Polvani, 2004. *An initial-value problem for testing numerical models of the global shallow-water equations*, Tellus A. DOI: [10.3402/tellusa.v56i5.14436](https://doi.org/10.3402/tellusa.v56i5.14436) -[^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. [10.1256/qj.06.12](https://doi.org/10.1256/qj.06.12) \ No newline at end of file +[^G04]: Galewsky, Scott, Polvani, 2004. *An initial-value problem for testing numerical models of the global shallow-water equations*, Tellus A. DOI: [10.3402/tellusa.v56i5.14436](https://doi.org/10.3402/tellusa.v56i5.14436) \ No newline at end of file diff --git a/docs/src/examples_3D.md b/docs/src/examples_3D.md new file mode 100644 index 000000000..d6b343b56 --- /dev/null +++ b/docs/src/examples_3D.md @@ -0,0 +1,160 @@ +# Examples 3D + +The following showcases several examples of SpeedyWeather.jl simulating +the [Primitive equations](@ref primitive_equation_model) with and without +humidity and with and without physical parameterizations. + +See also [Examples 2D](@ref Examples) for examples with the +[Barotropic vorticity equation](@ref barotropic_vorticity_model) and the +[shallow water model](@ref shallow_water_model). + +## Jablonowski-Williamson baroclinic wave + +```@example jablonowski +using SpeedyWeather +spectral_grid = SpectralGrid(trunc=31, nlev=8, Grid=FullGaussianGrid, dealiasing=3) + +orography = ZonalRidge(spectral_grid) +initial_conditions = InitialConditions( + vordiv = ZonalWind(), + temp = JablonowskiTemperature(), + pres = ZeroInitially()) + +model = PrimitiveDryModel(; spectral_grid, orography, initial_conditions, physics=false) +simulation = initialize!(model) +model.feedback.verbose = false # hide +run!(simulation, period=Day(9)) +nothing # hide +``` + +The Jablonowski-Williamson baroclinic wave test case[^JW06] using the +[Primitive equation model](@ref primitive_equation_model) particularly the dry model, +as we switch off all physics with `physics=false`. +We want to use 8 vertical levels, and a lower resolution of T31 on a +[full Gaussian grid](@ref FullGaussianGrid). +The Jablonowski-Williamson initial conditions are `ZonalWind` for vorticity and divergence +(curl and divergence of ``u, v``), `JablonowskiTemperature` for temperature and +`ZeroInitially` for pressure. The orography is just a `ZonalRidge`. +There is no forcing and the initial conditions are baroclinically unstable which kicks +off a wave propagating eastward. This wave becomes obvious when visualised with + +```@example jablonowski +using CairoMakie + +vor = simulation.diagnostic_variables.layers[end].grid_variables.vor_grid +heatmap(vor, title="Surface relative vorticity") +save("jablonowski.png", ans) # hide +nothing # hide +``` +![Jablonowski pyplot](jablonowski.png) + +## Held-Suarez forcing + + + +## Aquaplanet + +```@example aquaplanet +using SpeedyWeather + +# components +spectral_grid = SpectralGrid(trunc=31, nlev=5) +ocean = AquaPlanet(spectral_grid, temp_equator=302, temp_poles=273) +land_sea_mask = AquaPlanetMask(spectral_grid) +orography = NoOrography(spectral_grid) + +# create model, initialize, run +model = PrimitiveWetModel(; spectral_grid, ocean, land_sea_mask, orography) +simulation = initialize!(model) +model.feedback.verbose = false # hide +run!(simulation, period=Day(50)) +nothing # hide +``` + +Here we have defined an aquaplanet simulation by +- creating an `ocean::AquaPlanet`. This will use constant sea surface temperatures that only vary with latitude. +- creating a `land_sea_mask::AquaPlanetMask` this will use a land-sea mask with `false`=ocean everywhere. +- creating an `orography::NoOrography` which will have no orography and zero surface geopotential. + +All passed on to the model constructor for a `PrimitiveWetModel`, we have now a model with humidity +and physics parameterization as they are defined by default (typing `model` will give you an overview +of its components). We could have change the `model.land` and `model.vegetation` components too, +but given the land-sea masks masks those contributions to the surface fluxes anyway, this is not +necessary. Note that neither sea surface temperature, land-sea mask +or orography have to agree. It is possible to have an ocean on top of a mountain. +For an ocean grid-cell that is (partially) masked by the land-sea mask, its value will +be (fractionally) ignored in the calculation of surface fluxes (potentially leading +to a zero flux depending on land surface temperatures). + +Now with the following we visualize the surface humidity after the 50 days of +simulation. We use 50 days as without mountains it takes longer for the initial conditions to +become unstable. The surface humidity shows small-scale patches in the tropics, which is a result +of the convection scheme, causing updrafts and downdrafts in both humidity and temperature. + +```@example aquaplanet +using CairoMakie + +humid = simulation.diagnostic_variables.layers[end].grid_variables.humid_grid +heatmap(humid, title="Surface specific humidity [kg/kg]", colormap=:oslo) + +save("aquaplanet.png", ans) # hide +nothing # hide +``` +![Aquaplanet](aquaplanet.png) + +## Aquaplanet without (deep) convection + +Now we want to compare the previous simulation to a simulation without +deep convection, called `DryBettsMiller`, because it is the +[Betts-Miller convection](@ref BettsMiller) but with humidity set to zero +in which case the convection is always non-precipitating _shallow_ +(because the missing latent heat release from condensation makes it shallower) +convection. In fact, this convection is the default when using +the `PrimitiveDryModel`. Instead of redefining the [Aquaplanet](@ref) setup again, +we simply reuse these components `spectral_grid`, `ocean`, `land_sea_mask` +and `orography` (because `spectral_grid` hasn't changed this is possible). + +```@example aquaplanet +# Execute the code from Aquaplanet above first! +convection = DryBettsMiller(spectral_grid, time_scale=Hour(4)) + +# reuse other model components from before +model = PrimitiveWetModel(; spectral_grid, ocean, land_sea_mask, orography, convection) + +simulation = initialize!(model) +run!(simulation, period=Day(50)) + +humid = simulation.diagnostic_variables.layers[end].grid_variables.humid_grid +heatmap(humid, title="No deep convection: Surface specific humidity [kg/kg]", colormap=:oslo) +save("aquaplanet_nodeepconvection.png", ans) # hide +nothing # hide +``` + +But we also want to compare this to a setup where convection is completely +disabled, i.e. `convection = NoConvection()` (many of the `No` model components +don't require the `spectral_grid` to be passed on, but some do!) + +```@example aquaplanet +# Execute the code from Aquaplanet above first! +convection = NoConvection(spectral_grid) + +# reuse other model components from before +model = PrimitiveWetModel(; spectral_grid, ocean, land_sea_mask, orography, convection) + +simulation = initialize!(model) +run!(simulation, period=Day(50)) + +humid = simulation.diagnostic_variables.layers[end].grid_variables.humid_grid +heatmap(humid, title="No convection: Surface specific humidity [kg/kg]", colormap=:oslo) +save("aquaplanet_noconvection.png", ans) # hide +nothing # hide +``` + +And the comparison looks like + +![Aquaplanet, no deep convection](aquaplanet_nodeepconvection.png) +![Aquaplanet, no convection](aquaplanet_noconvection.png) + +## 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. [10.1256/qj.06.12](https://doi.org/10.1256/qj.06.12) \ No newline at end of file diff --git a/docs/src/forcing_drag.md b/docs/src/forcing_drag.md index aa51a1732..231de5b41 100644 --- a/docs/src/forcing_drag.md +++ b/docs/src/forcing_drag.md @@ -418,7 +418,7 @@ in grid-point space by writing into `u_tend_grid` and `v_tend_grid`. Now that we have defined a new forcing, as well as how to initialize it and what it is supposed to execute on every time step, we also want to use it. -We generally follow other [Examples](@ref), start with the [SpectralGrid](@ref) +We generally follow other [Examples](@ref Examples), start with the [SpectralGrid](@ref) and use that to get an instance of `StochasticStirring`. This calls the generator function from [Custom forcing: generator function](@ref). Here we want to stir vorticity not at the default latitude of 45N, but on the southern diff --git a/docs/src/how_to_run_speedy.md b/docs/src/how_to_run_speedy.md index 833b27811..07847f22b 100644 --- a/docs/src/how_to_run_speedy.md +++ b/docs/src/how_to_run_speedy.md @@ -11,7 +11,7 @@ user may want to adjust are chosen and live in their respective model components 4. [Run that simulation](@ref run). In the following we will describe these steps in more detail. -If you want to start with some examples first, see [Examples](@ref) +If you want to start with some examples first, see [Examples](@ref Examples) which has easy and more advanced examples of how to set up SpeedyWeather.jl to run the simulation you want. @@ -125,7 +125,7 @@ any name) ```@example howto model = ShallowWaterModel(; spectral_grid, time_stepping) ``` -This logic continues for all model components, see [Examples](@ref). +This logic continues for all model components, see [Examples](@ref Examples). All model components are also subtype (i.e. `<:`) of some abstract supertype, this feature can be made use of to redefine not just constant parameters of existing model components but also diff --git a/docs/src/installation.md b/docs/src/installation.md index 8ec0ecbfc..341e71b52 100644 --- a/docs/src/installation.md +++ b/docs/src/installation.md @@ -37,5 +37,5 @@ This is an [extension](https://pkgdocs.julialang.org/v1.10/creating-packages/#Conditional-loading-of-code-in-packages-(Extensions)), meaning that this functionality is only loaded from SpeedyWeather when `using Makie` (or its backends CairoMakie.jl, GLMakie.jl, ...). Hence, if you want to make use of this -extension (some [Examples](@ref) show this) you need to install Makie.jl manually. +extension (some [Examples](@ref Examples) show this) you need to install Makie.jl manually. From 68422b9b97ab043fde0b04d86f21202288634cef Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 21:19:36 -0400 Subject: [PATCH 02/14] NoTypes for surface fluxes --- src/models/primitive_dry.jl | 5 ++- src/models/primitive_wet.jl | 10 +++-- src/physics/surface_fluxes.jl | 83 ++++++++++++++++++++--------------- 3 files changed, 59 insertions(+), 39 deletions(-) diff --git a/src/models/primitive_dry.jl b/src/models/primitive_dry.jl index 6b2fb10b6..daae1491c 100644 --- a/src/models/primitive_dry.jl +++ b/src/models/primitive_dry.jl @@ -31,7 +31,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{ VD<:AbstractVerticalDiffusion, SUT<:AbstractSurfaceThermodynamics, SUW<:AbstractSurfaceWind, - SH<:AbstractSurfaceHeat, + SH<:AbstractSurfaceSensibleHeat, CV<:AbstractConvection, SW<:AbstractShortwave, LW<:AbstractLongwave, @@ -126,6 +126,9 @@ function initialize!(model::PrimitiveDry; time::DateTime = DEFAULT_DATE) initialize!(model.vertical_diffusion, model) initialize!(model.shortwave_radiation, model) initialize!(model.longwave_radiation, model) + initialize!(model.surface_thermodynamics, model) + initialize!(model.surface_wind, model) + initialize!(model.surface_heat_flux, model) # initial conditions prognostic_variables = PrognosticVariables(spectral_grid, model) diff --git a/src/models/primitive_wet.jl b/src/models/primitive_wet.jl index 802d590bd..595649721 100644 --- a/src/models/primitive_wet.jl +++ b/src/models/primitive_wet.jl @@ -34,8 +34,8 @@ Base.@kwdef mutable struct PrimitiveWetModel{ VD<:AbstractVerticalDiffusion, SUT<:AbstractSurfaceThermodynamics, SUW<:AbstractSurfaceWind, - SH<:AbstractSurfaceHeat, - EV<:AbstractEvaporation, + SH<:AbstractSurfaceSensibleHeat, + EV<:AbstractSurfaceEvaporation, LSC<:AbstractCondensation, CV<:AbstractConvection, SW<:AbstractShortwave, @@ -83,7 +83,7 @@ Base.@kwdef mutable struct PrimitiveWetModel{ surface_thermodynamics::SUT = SurfaceThermodynamicsConstant(spectral_grid) surface_wind::SUW = SurfaceWind(spectral_grid) surface_heat_flux::SH = SurfaceSensibleHeat(spectral_grid) - evaporation::EV = SurfaceEvaporation(spectral_grid) + surface_evaporation::EV = SurfaceEvaporation(spectral_grid) large_scale_condensation::LSC = ImplicitCondensation(spectral_grid) convection::CV = SimplifiedBettsMiller(spectral_grid) shortwave_radiation::SW = TransparentShortwave(spectral_grid) @@ -142,6 +142,10 @@ function initialize!(model::PrimitiveWet; time::DateTime = DEFAULT_DATE) initialize!(model.convection, model) initialize!(model.shortwave_radiation, model) initialize!(model.longwave_radiation, model) + initialize!(model.surface_thermodynamics, model) + initialize!(model.surface_wind, model) + initialize!(model.surface_heat_flux, model) + initialize!(model.surface_evaporation, model) # initial conditions prognostic_variables = PrognosticVariables(spectral_grid, model) diff --git a/src/physics/surface_fluxes.jl b/src/physics/surface_fluxes.jl index 92da45deb..b600341d6 100644 --- a/src/physics/surface_fluxes.jl +++ b/src/physics/surface_fluxes.jl @@ -1,28 +1,29 @@ -abstract type AbstractSurfaceWind <: AbstractParameterization end abstract type AbstractSurfaceThermodynamics <: AbstractParameterization end -abstract type AbstractSurfaceHeat <: AbstractParameterization end -abstract type AbstractEvaporation <: AbstractParameterization end +abstract type AbstractSurfaceWind <: AbstractParameterization end +abstract type AbstractSurfaceSensibleHeat <: AbstractParameterization end +abstract type AbstractSurfaceEvaporation <: AbstractParameterization end +# defines the order in which they are called und unpacks to dispatch function surface_fluxes!(column::ColumnVariables, model::PrimitiveEquation) # get temperature, humidity and density at surface - surface_thermodynamics!(column, model.surface_thermodynamics, model.atmosphere, model) + surface_thermodynamics!(column, model.surface_thermodynamics, model) # also calculates surface wind speed necessary for other fluxes too - surface_wind_stress!(column, model.surface_wind) + surface_wind_stress!(column, model.surface_wind, model) - # now call other heat and humidity fluxes - sensible_heat_flux!(column, model.surface_heat_flux, model.atmosphere) - evaporation!(column, model) + # now call other heat (wet and dry) and humidity fluxes (PrimitiveWet only) + sensible_heat_flux!(column, model.surface_heat_flux, model) + model isa PrimitiveWet && surface_evaporation!(column, model.surface_evaporation, model) end +## SURFACE THERMODYNAMICS export SurfaceThermodynamicsConstant -struct SurfaceThermodynamicsConstant{NF<:AbstractFloat} <: AbstractSurfaceThermodynamics end -SurfaceThermodynamicsConstant(SG::SpectralGrid) = SurfaceThermodynamicsConstant{SG.NF}() +struct SurfaceThermodynamicsConstant <: AbstractSurfaceThermodynamics end +SurfaceThermodynamicsConstant(SG::SpectralGrid) = SurfaceThermodynamicsConstant() function surface_thermodynamics!( column::ColumnVariables, ::SurfaceThermodynamicsConstant, - atmosphere::AbstractAtmosphere, model::PrimitiveWet) # surface value is same as lowest model level @@ -30,7 +31,7 @@ function surface_thermodynamics!( column::ColumnVariables, column.surface_humid = column.humid[end] # humidity at surface is the same as # surface air density via virtual temperature - (; R_dry) = atmosphere + (; R_dry) = model.atmosphere Tᵥ = column.temp_virt[column.nlev] column.surface_air_density = column.pres[end]/(R_dry*Tᵥ) end @@ -45,6 +46,13 @@ function surface_thermodynamics!( column::ColumnVariables, column.surface_air_density = column.pres[end]/(R_dry*column.surface_temp) end +## WIND STRESS +export NoSurfaceWind +struct NoSurfaceWind <: AbstractSurfaceWind end +NoSurfaceWind(::SpectralGrid) = NoSurfaceWind() +initialize!(::NoSurfaceWind, ::PrimitiveEquation) = nothing +surface_wind_stress!(::ColumnVariables, ::NoSurfaceWind, ::PrimitiveEquation) = nothing + export SurfaceWind Base.@kwdef struct SurfaceWind{NF<:AbstractFloat} <: AbstractSurfaceWind "Ratio of near-surface wind to lowest-level wind [1]" @@ -67,9 +75,11 @@ Base.@kwdef struct SurfaceWind{NF<:AbstractFloat} <: AbstractSurfaceWind end SurfaceWind(SG::SpectralGrid; kwargs...) = SurfaceWind{SG.NF}(; kwargs...) +initialize!(::SurfaceWind, ::PrimitiveEquation) = nothing -function surface_wind_stress!( column::ColumnVariables{NF}, - surface_wind::SurfaceWind) where NF +function surface_wind_stress!( column::ColumnVariables, + surface_wind::SurfaceWind, + model::PrimitiveEquation) (; land_fraction) = column (; f_wind, V_gust, drag_land, drag_sea) = surface_wind @@ -106,8 +116,15 @@ function surface_wind_stress!( column::ColumnVariables{NF}, return nothing end +## SENSIBLE HEAT FLUX +export NoSurfaceSensibleHeat +struct NoSurfaceSensibleHeat <: AbstractSurfaceSensibleHeat end +NoSurfaceSensibleHeat(::SpectralGrid) = NoSurfaceSensibleHeat() +initialize!(::NoSurfaceSensibleHeat, ::PrimitiveEquation) = nothing +sensible_heat_flux!(::ColumnVariables, ::NoSurfaceSensibleHeat, ::PrimitiveEquation) = nothing + export SurfaceSensibleHeat -Base.@kwdef struct SurfaceSensibleHeat{NF<:AbstractFloat} <: AbstractSurfaceHeat +Base.@kwdef struct SurfaceSensibleHeat{NF<:AbstractFloat} <: AbstractSurfaceSensibleHeat "Use (possibly) flow-dependent column.boundary_layer_drag coefficient" use_boundary_layer_drag::Bool = true @@ -123,13 +140,14 @@ Base.@kwdef struct SurfaceSensibleHeat{NF<:AbstractFloat} <: AbstractSurfaceHeat end SurfaceSensibleHeat(SG::SpectralGrid; kwargs...) = SurfaceSensibleHeat{SG.NF}(; kwargs...) +initialize!(::SurfaceSensibleHeat, ::PrimitiveEquation) = nothing function sensible_heat_flux!( column::ColumnVariables, heat_flux::SurfaceSensibleHeat, - atmosphere::AbstractAtmosphere + model::PrimitivEquation, ) - cₚ = atmosphere.heat_capacity + cₚ = model.atmosphere.heat_capacity (; heat_exchange_land, heat_exchange_sea, max_flux) = heat_flux ρ = column.surface_air_density @@ -173,12 +191,19 @@ function sensible_heat_flux!( return nothing end +## SURFACE EVAPORATION +export NoSurfaceEvaporation +struct NoSurfaceEvaporation <: AbstractSurfaceSensibleHeat end +NoSurfaceEvaporation(::SpectralGrid) = NoSurfaceEvaporation() +initialize!(::NoSurfaceEvaporation, ::PrimitiveEquation) = nothing +surface_evaporation!(::ColumnVariables, ::NoSurfaceEvaporation, ::PrimitiveEquation) = nothing + export SurfaceEvaporation """ Surface evaporation following a bulk formula with wind from model.surface_wind $(TYPEDFIELDS)""" -Base.@kwdef struct SurfaceEvaporation{NF<:AbstractFloat} <: AbstractEvaporation +Base.@kwdef struct SurfaceEvaporation{NF<:AbstractFloat} <: AbstractSurfaceEvaporation "Use column.boundary_layer_drag coefficient" use_boundary_layer_drag::Bool = true @@ -192,21 +217,9 @@ end SurfaceEvaporation(SG::SpectralGrid; kwargs...) = SurfaceEvaporation{SG.NF}(; kwargs...) -# don't do anything for dry core -function evaporation!( column::ColumnVariables, - model::PrimitiveDry) - return nothing -end - -# function barrier -function evaporation!( column::ColumnVariables, - model::PrimitiveWet) - evaporation!(column, model.evaporation, model.clausius_clapeyron) -end - -function evaporation!( column::ColumnVariables{NF}, - evaporation::SurfaceEvaporation, - clausius_clapeyron::AbstractClausiusClapeyron) where NF +function surface_evaporation!( column::ColumnVariables, + evaporation::SurfaceEvaporation, + model::PrimitiveWet) (; skin_temperature_sea, skin_temperature_land, pres) = column (; moisture_exchange_land, moisture_exchange_sea) = evaporation @@ -214,8 +227,8 @@ function evaporation!( column::ColumnVariables{NF}, # SATURATION HUMIDITY OVER LAND AND OCEAN surface_pressure = pres[end] - sat_humid_land = saturation_humidity(skin_temperature_land, surface_pressure, clausius_clapeyron) - sat_humid_sea = saturation_humidity(skin_temperature_sea, surface_pressure, clausius_clapeyron) + sat_humid_land = saturation_humidity(skin_temperature_land, surface_pressure, model.clausius_clapeyron) + sat_humid_sea = saturation_humidity(skin_temperature_sea, surface_pressure, model.clausius_clapeyron) ρ = column.surface_air_density V₀ = column.surface_wind_speed From ebcaea425b273d7def9edd3cdce0a598f42185dd Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 21:19:57 -0400 Subject: [PATCH 03/14] Held-Suarez example added --- README.md | 2 +- docs/src/convection.md | 3 ++- docs/src/examples_3D.md | 59 ++++++++++++++++++++++++++++++++++++++++- 3 files changed, 61 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index d52d62953..ac9912bcb 100644 --- a/README.md +++ b/README.md @@ -73,7 +73,7 @@ about dos and don'ts. Just express your interest to contribute and we'll be happ ## Example use For a more comprehensive tutorial with several examples, see -[Examples](https://speedyweather.github.io/SpeedyWeather.jl/dev/examples/) in the documentation. +[Examples](https://speedyweather.github.io/SpeedyWeather.jl/dev/examples_2D/) in the documentation. The interface to SpeedyWeather.jl consist of 5 steps: define the grid, create model components, construct the model, initialize, run diff --git a/docs/src/convection.md b/docs/src/convection.md index eb452fbf7..e1a87ea74 100644 --- a/docs/src/convection.md +++ b/docs/src/convection.md @@ -152,7 +152,8 @@ Q_{ref} &= \int_{p_0}^{p_{LZB}} -q_{ref} dp \\ f_q &= 1 - \frac{\Delta q}{Q_ref} \\ q_{ref, 2} &= f_q q_{ref} \\ \Delta T &= \frac{1}{\Delta p} \int_{p_0}^{p_{LZB}} -(T - T_{ref}) dp \\ -T_{ref,2} = T_{ref} - \Delta T +T_{ref,2} &= T_{ref} - \Delta T +\end{aligned} ``` ## Corrected relaxation diff --git a/docs/src/examples_3D.md b/docs/src/examples_3D.md index d6b343b56..0efaa88bc 100644 --- a/docs/src/examples_3D.md +++ b/docs/src/examples_3D.md @@ -50,7 +50,62 @@ nothing # hide ## Held-Suarez forcing +```@example heldsuarez +using SpeedyWeather +spectral_grid = SpectralGrid(trunc=31, nlev=8) + +# construct model with only Held-Suarez forcing, no other physics +model = PrimitiveDryModel(; + spectral_grid, + + # Held-Suarez forcing and drag + temperature_relaxation = HeldSuarez(spectral_grid), + drag = LinearDrag(spectral_grid), + + # switch off other physics + convection = NoConvection(), + shortwave_radiation = NoShortwave(), + longwave_radiation = NoLongwave(), + vertical_diffusion = NoVerticalDiffusion(), + + # switch off surface fluxes (makes ocean/land/land-sea mask redundant) + surface_wind = NoSurfaceWind(), + surface_evaporation = NoSurfaceEvaporation(), + surface_sensible_heat = NoSurfaceSensibleHeat(), + # use Earth's orography + orography = EarthOrography(spectral_grid) +) + +simulation = initialize!(model) +model.feedback.verbose = false # hide +run!(simulation, period=Day(20)) +nothing # hide +``` + +The code above defines the Held-Suarez forcing [^HS94] in terms of temperature relaxation +and a linear drag term that is applied near the planetary boundary but switches off +all other physics in the primitive equation model without humidity. +Switching off the surface wind would also automatically turn off the surface evaporation +and sensible heat flux as that one is proportional to the surface wind (which is zero +with `NoSurfaceWind`). But to also avoid the calculation being run at all we use +`NoSurfaceEvaporation()` and `NoSurfaceSensibleHeat()` for the model constructor. +Many of the `NoSomething` model components do not require the +spectral grid to be passed on, but as a convention we allow every model component +to have it for construction even if not required. + +Visualising surface temperature with + +```@example heldsuarez +using CairoMakie + +temp = simulation.diagnostic_variables.layers[end].grid_variables.temp_grid +heatmap(temp, title="Surface temperature [K]", colormap=:thermal) + +save("heldsuarez.png", ans) # hide +nothing # hide +``` +![Held-Suarez](heldsuarez.png) ## Aquaplanet @@ -157,4 +212,6 @@ And the comparison looks like ## 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. [10.1256/qj.06.12](https://doi.org/10.1256/qj.06.12) \ No newline at end of file +[^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) + +[^HS94]: Held, I. M. & Suarez, M. J. A Proposal for the Intercomparison of the Dynamical Cores of Atmospheric General Circulation Models. Bulletin of the American Meteorological Society 75, 1825-1830 (1994). DOI:[10.1175/1520-0477(1994)075<1825:APFTIO>2.0.CO;2](https://doi.org/10.1175/1520-0477(1994)075<1825:APFTIO>2.0.CO;2) From 9f44324fbce79e348dc13a41c8048024b96ceb59 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 21:23:36 -0400 Subject: [PATCH 04/14] Gremlins removed in temperature_relaxation.jl --- src/physics/temperature_relaxation.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/physics/temperature_relaxation.jl b/src/physics/temperature_relaxation.jl index f9f37df33..b3eb6443f 100644 --- a/src/physics/temperature_relaxation.jl +++ b/src/physics/temperature_relaxation.jl @@ -14,7 +14,7 @@ temperature_relaxation!(::ColumnVariables, ::NoTemperatureRelaxation, ::Primitiv export HeldSuarez """ -Struct that defines the temperature relaxation from Held and Suarez, 1996 BAMS +Temperature relaxation from Held and Suarez, 1996 BAMS $(TYPEDFIELDS)""" Base.@kwdef struct HeldSuarez{NF<:AbstractFloat} <: AbstractTemperatureRelaxation # DIMENSIONS @@ -50,9 +50,9 @@ Base.@kwdef struct HeldSuarez{NF<:AbstractFloat} <: AbstractTemperatureRelaxatio κ::Base.RefValue{NF} = Ref(zero(NF)) p₀::Base.RefValue{NF} = Ref(zero(NF)) - temp_relax_freq::Matrix{NF} = zeros(NF, nlev, nlat) # (inverse) relax time scale per layer and lat - temp_equil_a::Vector{NF} = zeros(NF, nlat) # terms to calc equilibrium temper func - temp_equil_b::Vector{NF} = zeros(NF, nlat) # of latitude and pressure + temp_relax_freq::Matrix{NF} = zeros(NF, nlev, nlat) # (inverse) relax time scale per layer and lat + temp_equil_a::Vector{NF} = zeros(NF, nlat) # terms to calc equilibrium temper func + temp_equil_b::Vector{NF} = zeros(NF, nlat) # of latitude and pressure end """ @@ -111,7 +111,7 @@ function temperature_relaxation!( atmosphere::AbstractAtmosphere, ) (; temp, temp_tend, pres, ln_pres) = column - j = column.jring[] # latitude ring index j + j = column.jring[] # latitude ring index j (; Tmin, temp_relax_freq, temp_equil_a, temp_equil_b) = scheme # surface reference pressure [Pa] and thermodynamic kappa R_dry/cₚ @@ -165,7 +165,7 @@ Base.@kwdef mutable struct JablonowskiRelaxation{NF<:AbstractFloat} <: AbstractT # precomputed constants, allocate here, fill in initialize! temp_relax_freq::Matrix{NF} = zeros(NF, nlev, nlat) # (inverse) relax time scale per layer and lat - temp_equil::Matrix{NF} = zeros(NF, nlev, nlat) # terms to calc equilibrium temperature as func + temp_equil::Matrix{NF} = zeros(NF, nlev, nlat) # terms to calc equilibrium temperature as func end """ @@ -236,7 +236,7 @@ function temperature_relaxation!( column::ColumnVariables, scheme::JablonowskiRelaxation) (; temp, temp_tend) = column - j = column.jring[] # latitude ring index j + j = column.jring[] # latitude ring index j (; temp_relax_freq, temp_equil) = scheme @inbounds for k in eachlayer(column) @@ -244,6 +244,6 @@ function temperature_relaxation!( column::ColumnVariables, # Held and Suarez 1996, equation 2, but using temp_equil from # Jablonowski and Williamson 2006, equation 6 - temp_tend[k] -= kₜ*(temp[k] - temp_equil[k, j]) + temp_tend[k] -= kₜ*(temp[k] - temp_equil[k, j]) end end From 8c486c5f60249561edf0144d08d0997f32eb2fc0 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 21:27:01 -0400 Subject: [PATCH 05/14] hide InteractiveUtils --- docs/src/radiation.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/radiation.md b/docs/src/radiation.md index e46342dfa..5902e1d7e 100644 --- a/docs/src/radiation.md +++ b/docs/src/radiation.md @@ -5,7 +5,8 @@ Currently implemented is ```@example radiation -using SpeedyWeather, InteractiveUtils +using InteractiveUtils # hide +using SpeedyWeather subtypes(SpeedyWeather.AbstractLongwave) ``` From 656f5fe05ada2aceca711a20e15325964fab0db9 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 21:46:17 -0400 Subject: [PATCH 06/14] Rename to SurfaceHeatFlux, surface_heat_flux, no sensible anymore --- docs/src/examples_3D.md | 12 ++++++------ src/models/primitive_dry.jl | 4 ++-- src/models/primitive_wet.jl | 4 ++-- src/physics/surface_fluxes.jl | 32 ++++++++++++++++---------------- 4 files changed, 26 insertions(+), 26 deletions(-) diff --git a/docs/src/examples_3D.md b/docs/src/examples_3D.md index 0efaa88bc..8bfe94c3f 100644 --- a/docs/src/examples_3D.md +++ b/docs/src/examples_3D.md @@ -60,7 +60,7 @@ model = PrimitiveDryModel(; # Held-Suarez forcing and drag temperature_relaxation = HeldSuarez(spectral_grid), - drag = LinearDrag(spectral_grid), + boundary_layer_drag = LinearDrag(spectral_grid), # switch off other physics convection = NoConvection(), @@ -70,8 +70,7 @@ model = PrimitiveDryModel(; # switch off surface fluxes (makes ocean/land/land-sea mask redundant) surface_wind = NoSurfaceWind(), - surface_evaporation = NoSurfaceEvaporation(), - surface_sensible_heat = NoSurfaceSensibleHeat(), + surface_heat_flux = NoSurfaceHeatFlux(), # use Earth's orography orography = EarthOrography(spectral_grid) @@ -87,9 +86,10 @@ The code above defines the Held-Suarez forcing [^HS94] in terms of temperature r and a linear drag term that is applied near the planetary boundary but switches off all other physics in the primitive equation model without humidity. Switching off the surface wind would also automatically turn off the surface evaporation -and sensible heat flux as that one is proportional to the surface wind (which is zero -with `NoSurfaceWind`). But to also avoid the calculation being run at all we use -`NoSurfaceEvaporation()` and `NoSurfaceSensibleHeat()` for the model constructor. +(not relevant in the primitive _dry_ model) and sensible heat flux as that one is proportional +to the surface wind (which is zero with `NoSurfaceWind`). +But to also avoid the calculation being run at all we use `NoSurfaceHeatFlux()` +for the model constructor. Many of the `NoSomething` model components do not require the spectral grid to be passed on, but as a convention we allow every model component to have it for construction even if not required. diff --git a/src/models/primitive_dry.jl b/src/models/primitive_dry.jl index daae1491c..01d710d33 100644 --- a/src/models/primitive_dry.jl +++ b/src/models/primitive_dry.jl @@ -31,7 +31,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{ VD<:AbstractVerticalDiffusion, SUT<:AbstractSurfaceThermodynamics, SUW<:AbstractSurfaceWind, - SH<:AbstractSurfaceSensibleHeat, + SH<:AbstractSurfaceHeatFlux, CV<:AbstractConvection, SW<:AbstractShortwave, LW<:AbstractLongwave, @@ -73,7 +73,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{ vertical_diffusion::VD = BulkRichardsonDiffusion(spectral_grid) surface_thermodynamics::SUT = SurfaceThermodynamicsConstant(spectral_grid) surface_wind::SUW = SurfaceWind(spectral_grid) - surface_heat_flux::SH = SurfaceSensibleHeat(spectral_grid) + surface_heat_flux::SH = SurfaceHeatFlux(spectral_grid) convection::CV = DryBettsMiller(spectral_grid) shortwave_radiation::SW = TransparentShortwave(spectral_grid) longwave_radiation::LW = JeevanjeeRadiation(spectral_grid) diff --git a/src/models/primitive_wet.jl b/src/models/primitive_wet.jl index 595649721..7a4c5eda8 100644 --- a/src/models/primitive_wet.jl +++ b/src/models/primitive_wet.jl @@ -34,7 +34,7 @@ Base.@kwdef mutable struct PrimitiveWetModel{ VD<:AbstractVerticalDiffusion, SUT<:AbstractSurfaceThermodynamics, SUW<:AbstractSurfaceWind, - SH<:AbstractSurfaceSensibleHeat, + SH<:AbstractSurfaceHeatFlux, EV<:AbstractSurfaceEvaporation, LSC<:AbstractCondensation, CV<:AbstractConvection, @@ -82,7 +82,7 @@ Base.@kwdef mutable struct PrimitiveWetModel{ vertical_diffusion::VD = BulkRichardsonDiffusion(spectral_grid) surface_thermodynamics::SUT = SurfaceThermodynamicsConstant(spectral_grid) surface_wind::SUW = SurfaceWind(spectral_grid) - surface_heat_flux::SH = SurfaceSensibleHeat(spectral_grid) + surface_heat_flux::SH = SurfaceHeatFlux(spectral_grid) surface_evaporation::EV = SurfaceEvaporation(spectral_grid) large_scale_condensation::LSC = ImplicitCondensation(spectral_grid) convection::CV = SimplifiedBettsMiller(spectral_grid) diff --git a/src/physics/surface_fluxes.jl b/src/physics/surface_fluxes.jl index b600341d6..40c1f0cc6 100644 --- a/src/physics/surface_fluxes.jl +++ b/src/physics/surface_fluxes.jl @@ -1,6 +1,6 @@ abstract type AbstractSurfaceThermodynamics <: AbstractParameterization end abstract type AbstractSurfaceWind <: AbstractParameterization end -abstract type AbstractSurfaceSensibleHeat <: AbstractParameterization end +abstract type AbstractSurfaceHeatFlux <: AbstractParameterization end abstract type AbstractSurfaceEvaporation <: AbstractParameterization end # defines the order in which they are called und unpacks to dispatch @@ -13,7 +13,7 @@ function surface_fluxes!(column::ColumnVariables, model::PrimitiveEquation) surface_wind_stress!(column, model.surface_wind, model) # now call other heat (wet and dry) and humidity fluxes (PrimitiveWet only) - sensible_heat_flux!(column, model.surface_heat_flux, model) + surface_heat_flux!(column, model.surface_heat_flux, model) model isa PrimitiveWet && surface_evaporation!(column, model.surface_evaporation, model) end @@ -117,14 +117,14 @@ function surface_wind_stress!( column::ColumnVariables, end ## SENSIBLE HEAT FLUX -export NoSurfaceSensibleHeat -struct NoSurfaceSensibleHeat <: AbstractSurfaceSensibleHeat end -NoSurfaceSensibleHeat(::SpectralGrid) = NoSurfaceSensibleHeat() -initialize!(::NoSurfaceSensibleHeat, ::PrimitiveEquation) = nothing -sensible_heat_flux!(::ColumnVariables, ::NoSurfaceSensibleHeat, ::PrimitiveEquation) = nothing - -export SurfaceSensibleHeat -Base.@kwdef struct SurfaceSensibleHeat{NF<:AbstractFloat} <: AbstractSurfaceSensibleHeat +export NoSurfaceHeatFlux +struct NoSurfaceHeatFlux <: AbstractSurfaceHeatFlux end +NoSurfaceHeatFlux(::SpectralGrid) = NoSurfaceHeatFlux() +initialize!(::NoSurfaceHeatFlux, ::PrimitiveEquation) = nothing +surface_heat_flux!(::ColumnVariables, ::NoSurfaceHeatFlux, ::PrimitiveEquation) = nothing + +export SurfaceHeatFlux +Base.@kwdef struct SurfaceHeatFlux{NF<:AbstractFloat} <: AbstractSurfaceHeatFlux "Use (possibly) flow-dependent column.boundary_layer_drag coefficient" use_boundary_layer_drag::Bool = true @@ -139,13 +139,13 @@ Base.@kwdef struct SurfaceSensibleHeat{NF<:AbstractFloat} <: AbstractSurfaceSens max_flux::NF = 100 end -SurfaceSensibleHeat(SG::SpectralGrid; kwargs...) = SurfaceSensibleHeat{SG.NF}(; kwargs...) -initialize!(::SurfaceSensibleHeat, ::PrimitiveEquation) = nothing +SurfaceHeatFlux(SG::SpectralGrid; kwargs...) = SurfaceHeatFlux{SG.NF}(; kwargs...) +initialize!(::SurfaceHeatFlux, ::PrimitiveEquation) = nothing -function sensible_heat_flux!( +function surface_heat_flux!( column::ColumnVariables, - heat_flux::SurfaceSensibleHeat, - model::PrimitivEquation, + heat_flux::SurfaceHeatFlux, + model::PrimitiveEquation, ) cₚ = model.atmosphere.heat_capacity (; heat_exchange_land, heat_exchange_sea, max_flux) = heat_flux @@ -193,7 +193,7 @@ end ## SURFACE EVAPORATION export NoSurfaceEvaporation -struct NoSurfaceEvaporation <: AbstractSurfaceSensibleHeat end +struct NoSurfaceEvaporation <: AbstractSurfaceHeatFlux end NoSurfaceEvaporation(::SpectralGrid) = NoSurfaceEvaporation() initialize!(::NoSurfaceEvaporation, ::PrimitiveEquation) = nothing surface_evaporation!(::ColumnVariables, ::NoSurfaceEvaporation, ::PrimitiveEquation) = nothing From 10e2c84ddfa137025eece98a67048d3a96c7d538 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 21:55:51 -0400 Subject: [PATCH 07/14] initialize! function was missing --- src/physics/surface_fluxes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/surface_fluxes.jl b/src/physics/surface_fluxes.jl index 40c1f0cc6..fb2577250 100644 --- a/src/physics/surface_fluxes.jl +++ b/src/physics/surface_fluxes.jl @@ -21,6 +21,7 @@ end export SurfaceThermodynamicsConstant struct SurfaceThermodynamicsConstant <: AbstractSurfaceThermodynamics end SurfaceThermodynamicsConstant(SG::SpectralGrid) = SurfaceThermodynamicsConstant() +initialize!(::SurfaceThermodynamicsConstant,::PrimitiveEquation) = nothing function surface_thermodynamics!( column::ColumnVariables, ::SurfaceThermodynamicsConstant, @@ -38,9 +39,8 @@ end function surface_thermodynamics!( column::ColumnVariables, ::SurfaceThermodynamicsConstant, - atmosphere::AbstractAtmosphere, model::PrimitiveDry) - (; R_dry) = atmosphere + (; R_dry) = model.atmosphere # surface value is same as lowest model level column.surface_temp = column.temp[end] # todo use constant POTENTIAL temperature column.surface_air_density = column.pres[end]/(R_dry*column.surface_temp) From ad2f017474d3aa714789f891231e861f08348885 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 22:01:53 -0400 Subject: [PATCH 08/14] initialize! for surface evaporation wasn't defined --- src/physics/surface_fluxes.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/physics/surface_fluxes.jl b/src/physics/surface_fluxes.jl index fb2577250..d8f78c36e 100644 --- a/src/physics/surface_fluxes.jl +++ b/src/physics/surface_fluxes.jl @@ -216,6 +216,7 @@ Base.@kwdef struct SurfaceEvaporation{NF<:AbstractFloat} <: AbstractSurfaceEvapo end SurfaceEvaporation(SG::SpectralGrid; kwargs...) = SurfaceEvaporation{SG.NF}(; kwargs...) +initialize!(::SurfaceEvaporation, ::PrimitiveWet) = nothing function surface_evaporation!( column::ColumnVariables, evaporation::SurfaceEvaporation, From b527135d74e00cef44fbeea7bd09c8a9098d3554 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 23:06:37 -0400 Subject: [PATCH 09/14] doc link corrected --- docs/src/convection.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/convection.md b/docs/src/convection.md index e1a87ea74..3e07d35f1 100644 --- a/docs/src/convection.md +++ b/docs/src/convection.md @@ -15,7 +15,7 @@ precipitating as the condensed humidity forms cloud droplets that eventually fall down as convective precipitation. See also [Large-scale condensation](@ref) in comparison. -## [Simplified Betts-Miller convective adjustment](@ref BettsMiller) +## [Simplified Betts-Miller convective adjustment](@id BettsMiller) We follow the simplification of the Betts-Miller convection scheme [^Betts1986][^BettsMiller1986] as studied by Frierson, 2007 [^Frierson2007]. From 1ff121508b6919b0cba0c7ec43f8cf552df71f74 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 23:22:36 -0400 Subject: [PATCH 10/14] Current implementations for condensation/convection automatically listed --- docs/src/convection.md | 26 +++++++++++++++++++++++--- docs/src/large_scale_condensation.md | 26 ++++++++++++++++++++++---- 2 files changed, 45 insertions(+), 7 deletions(-) diff --git a/docs/src/convection.md b/docs/src/convection.md index 3e07d35f1..eba7690d2 100644 --- a/docs/src/convection.md +++ b/docs/src/convection.md @@ -15,7 +15,17 @@ precipitating as the condensed humidity forms cloud droplets that eventually fall down as convective precipitation. See also [Large-scale condensation](@ref) in comparison. -## [Simplified Betts-Miller convective adjustment](@id BettsMiller) +## Convection implementations + +Currently implemented are +```@example convection +using InteractiveUtils # hide +using SpeedyWeather +subtypes(SpeedyWeather.AbstractConvection) +``` +which are described in the following. + +## [Simplified Betts-Miller convection](@id BettsMiller) We follow the simplification of the Betts-Miller convection scheme [^Betts1986][^BettsMiller1986] as studied by Frierson, 2007 [^Frierson2007]. @@ -24,9 +34,9 @@ convection as an adjustment towards a (pseudo-) moist adiabat reference profile and its associated humidity profile. Meaning that conceptually for every vertical column in the atmosphere we -1. Diagnose the vertical temperature and humidity profile of the environment relative to the surface and calculate the adiabat to the level of zero buoyancy. +1. Diagnose the vertical temperature and humidity profile of the environment relative to the adiabat up to the level of zero buoyancy. 2. Decide whether convection should take place and whether it is deep (precipitating) or shallow (non-precipitating). -3. Relax temperature and humidity towards (adjusted) profiles from 1. +3. Relax temperature and humidity towards (corrected) profiles from 1. ## Reference profiles @@ -183,6 +193,16 @@ P = -\int \frac{\Delta t}{g \rho} \delta q dp In the shallow convection case ``P=0`` due to the correction even though in the first guess relaxation ``P<0`` was possible, but for deep convection ``P>0`` by definition. +## Dry convection + +In the primitive equation model with humidity the [Betts-Miller convection scheme](@ref BettsMiller) +as described above is defined. Without humidity, a dry version reduces to the +[Shallow convection](@ref) case. The two different shallow convection schemes in +Frierson 2007[^Frierson2007], the "shallower" shallow convection scheme and the "qref" +(as implemented here in [Shallow convection](@ref)) in that case also reduce to +the same formulation. The dry Betts-Miller convection scheme is the default +in the primitive equation model without humidity. + ## References [^Betts1986]: Betts, A. K., 1986: A new convective adjustment scheme. Part I: Observational and theoretical basis. Quart. J. Roy. Meteor. Soc.,112, 677-691. DOI: [10.1002/qj.49711247307](https://doi.org/10.1002/qj.49711247307) diff --git a/docs/src/large_scale_condensation.md b/docs/src/large_scale_condensation.md index 3a3138ff8..8a1dd3290 100644 --- a/docs/src/large_scale_condensation.md +++ b/docs/src/large_scale_condensation.md @@ -11,6 +11,20 @@ quantities such as specific humidity, pressure and temperature within a given grid cell, even though there might be considerably variability of these quantities within a grid cell if the resolution was higher. +## Condensation implementations + +Currently implemented are + +```@example condensation +using InteractiveUtils # hide +using SpeedyWeather +subtypes(SpeedyWeather.AbstractCondensation) +``` + +which are described in the following. + +## Explicit large-scale condensation + We parameterize this process of _large-scale condensation_ when relative humidity in a grid cell reaches saturation and remove the excess humidity quickly (given time integration constraints, see below) and with an implicit (in the time @@ -32,10 +46,14 @@ quantities at time step ``i`` to calculate the tendency. The latent heat release condensation is in the second equation. However, treating this explicitly poses the problem that because the saturation humidity is calculated from the current temperature ``T_i``, which is increased due to the latent heat release, the humidity after this time step will be -undersaturated. Ideally, one would want to condense towards the new saturation humidity -``q^\star(T_{i+1})`` so that condensation draws the relative humidity back down to 100% not below it. -Taylor expansion at ``i`` with ``\Delta T = T_{i+1} - T_i`` (and ``\delta q`` similarly) -to first order yields +undersaturated. + +## Implicit large-scale condensation + +Ideally, one would want to condense towards the _new_ saturation humidity +``q^\star(T_{i+1})`` at ``i+1`` so that condensation draws the relative humidity back down +to 100% not below it. Taylor expansion at ``i`` of the equation above with ``q^\star(T_{i+1})`` +and ``\Delta T = T_{i+1} - T_i`` (and ``\Delta q`` similarly) to first order yields ```math q_{i+1} - q_i = q^\star(T_{i+1}) - q_i = q^\star(T_i) + (T_{i+1} - T_i) From 7d9d6330f91005c98d4ccf59f4cf4999f61b5a93 Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 4 Apr 2024 23:35:55 -0400 Subject: [PATCH 11/14] Last tweaks --- docs/src/examples_2D.md | 19 +++++++++++-------- docs/src/examples_3D.md | 2 ++ docs/src/output.md | 2 +- 3 files changed, 14 insertions(+), 9 deletions(-) diff --git a/docs/src/examples_2D.md b/docs/src/examples_2D.md index 3b968dd0d..601b22293 100644 --- a/docs/src/examples_2D.md +++ b/docs/src/examples_2D.md @@ -139,9 +139,12 @@ nothing # hide You see that in comparison the unicode plot heavily coarse-grains the simulation, well it's unicode after all! Here, we have unpacked the netCDF file using [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl) and then plotted via `heatmap(lon, lat, vor)`. While you can do that to give you more control -on the plotting, SpeedyWeather.jl also defines an extension for Makie.jl, see [Extensions](@ref) -because if our matrix `vor` here was an `AbstractGrid` (see [RingGrids](@ref)) then all -its geographic information (which grid point is where) is encoded in the type. So we can also do +on the plotting, SpeedyWeather.jl also defines an extension for Makie.jl, see [Extensions](@ref). +Because if our matrix `vor` here was an `AbstractGrid` (see [RingGrids](@ref)) then all +its geographic information (which grid point is where) would directly be encoded in the type. +From the netCDF file you need to use the longitude and latitude dimensions. + +So we can also just do ```@example galewsky_setup vor_grid = FullGaussianGrid(vor) @@ -154,7 +157,8 @@ nothing # hide Note that here you need to know which grid the data comes on (an error is thrown if `FullGaussianGrid(vor)` is not size compatible). By default the output will be on the FullGaussianGrid, but if you -play around with other grids, you'd need to change this here. +play around with other grids, you'd need to change this here, +see [NetCDF output on other grids](@ref output_grid). We did want to showcase the usage of [NetCDF output](@ref) here, but from now on we will use `heatmap` to plot data on our grids directly, without storing output first. @@ -162,7 +166,6 @@ So for our current simulation, that means at time = 12 days, vorticity on the gr is stored in the diagnostic variables and can be visualised with ```@example galewsky_setup -t = ds.dim["time"] vor = simulation.diagnostic_variables.layers[1].grid_variables.vor_grid heatmap(vor, title="Relative vorticity [1/s]") save("galewsky2.png", ans) # hide @@ -178,9 +181,9 @@ Let's try it out! We create an `EarthOrography` struct like so orography = EarthOrography(spectral_grid) ``` -It will read the orography from file as shown, and there are some smoothing options too, but let's not change them. -Same as before, create a model, initialize into a simulation, run. This time directly for 12 days so that we can -compare with the last plot +It will read the orography from file as shown (only at `initialize!(model)`), and there are some smoothing +options too, but let's not change them. Same as before, create a model, initialize into a simulation, run. +This time directly for 12 days so that we can compare with the last plot ```@example galewsky_setup model = ShallowWaterModel(; spectral_grid, orography, initial_conditions) diff --git a/docs/src/examples_3D.md b/docs/src/examples_3D.md index 8bfe94c3f..1f0293c8c 100644 --- a/docs/src/examples_3D.md +++ b/docs/src/examples_3D.md @@ -177,6 +177,7 @@ convection = DryBettsMiller(spectral_grid, time_scale=Hour(4)) model = PrimitiveWetModel(; spectral_grid, ocean, land_sea_mask, orography, convection) simulation = initialize!(model) +model.feedback.verbose = false # hide run!(simulation, period=Day(50)) humid = simulation.diagnostic_variables.layers[end].grid_variables.humid_grid @@ -197,6 +198,7 @@ convection = NoConvection(spectral_grid) model = PrimitiveWetModel(; spectral_grid, ocean, land_sea_mask, orography, convection) simulation = initialize!(model) +model.feedback.verbose = false # hide run!(simulation, period=Day(50)) humid = simulation.diagnostic_variables.layers[end].grid_variables.humid_grid diff --git a/docs/src/output.md b/docs/src/output.md index 91e2798c3..6c72287a4 100644 --- a/docs/src/output.md +++ b/docs/src/output.md @@ -86,7 +86,7 @@ ds["time"][:] ``` very neatly hourly output in the NetCDF file! -## Example 2: Output onto a higher/lower resolution grid +## [Example 2: Output onto a higher/lower resolution grid](@id output_grid) Say we want to run the model at a given horizontal resolution but want to output on another resolution, the `OutputWriter` takes as argument `output_Grid<:AbstractFullGrid` and `nlat_half::Int`. From 568af0fa3101592cea776c40d5b608fdcad5157d Mon Sep 17 00:00:00 2001 From: Milan Date: Fri, 5 Apr 2024 00:28:07 -0400 Subject: [PATCH 12/14] SpeedyTransforms docs with Makie too --- docs/src/speedytransforms.md | 44 +++++++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 11 deletions(-) diff --git a/docs/src/speedytransforms.md b/docs/src/speedytransforms.md index 71b839336..730ed6662 100644 --- a/docs/src/speedytransforms.md +++ b/docs/src/speedytransforms.md @@ -186,8 +186,14 @@ We now wrap this matrix therefore to associate it with the necessary grid information ```@example speedytransforms map = FullClenshawGrid(m) -plot(map) + +using CairoMakie +heatmap(map) +save("random_pattern.png", ans) # hide +nothing # hide ``` +![Random pattern](random_pattern.png) + Now we transform into spectral space and call `power_spectrum(::LowerTriangularMatrix)` ```@example speedytransforms alms = spectral(map) @@ -225,12 +231,17 @@ k = 1:32 alms = randn(LowerTriangularMatrix{Complex{Float32}}, 32, 32) alms .*= k.^-2 ``` -Awesome. For higher degrees and orders the amplitude clearly decreases! Now -to grid-point space and let us visualize the result +Awesome. For higher degrees and orders the amplitude clearly decreases! +Now to grid-point space and let us visualize the result ```@example speedytransforms map = gridded(alms) -plot(map) + +using CairoMakie +heatmap(map, title="k⁻²-distributed noise") +save("random_noise.png", ans) # hide +nothing # hide ``` +!(Random noise)[random_noise.png] You can always access the underlying data in `map` via `map.data` in case you need to get rid of the wrapping into a grid again! @@ -348,7 +359,7 @@ S = SpectralTransform(u, one_more_degree=true) us = spectral(u, S) vs = spectral(v, S) -vor = curl(us, vs) +vor = curl(us, vs) / spectral_grid.radius ``` (Copies of) the velocity fields are unscaled by the cosine of latitude (see above), then transformed into spectral space, and the returned `vor` requires a manual division @@ -388,20 +399,31 @@ Now we need to apply the inverse Laplace operator to ``f\zeta/g`` which we do as ```@example speedytransforms fζ_g_spectral = spectral(fζ_g, one_more_degree=true) -η = SpeedyTransforms.∇⁻²(fζ_g_spectral) * spectral_grid.radius^2 + +R = spectral_grid.radius +η = SpeedyTransforms.∇⁻²(fζ_g_spectral) * R^2 η_grid = gridded(η, Grid=spectral_grid.Grid) nothing # hide ``` - Note the manual scaling with the radius ``R^2`` here. We now compare the results ```@example speedytransforms -plot(η_grid) +using CairoMakie +heatmap(η_grid, title="Geostrophic interface displacement η [m]") +save("eta_geostrophic.png", ans) # hide +nothing # hide ``` -Which is the interface displacement assuming geostrophy. The actual interface -displacement contains also ageostrophy +!(Geostrophic eta)[eta_geostrophic.png] + +Which is the interface displacement assuming geostrophy. +The actual interface displacement contains also ageostrophy ```@example speedytransforms -plot(simulation.diagnostic_variables.surface.pres_grid) +η_grid2 = simulation.diagnostic_variables.surface.pres_grid +heatmap(η_grid2, title="Interface displacement η [m] with ageostrophy") +save("eta_ageostrophic.png", ans) # hide +nothing # hide ``` +!(Ageostrophic eta)[eta_ageostrophic.png] + Strikingly similar! The remaining differences are the ageostrophic motions but also note that the mean is off. This is because geostrophy only use/defines the gradient of ``\eta`` not the absolute values itself. Our geostrophic ``\eta_g`` has by construction From d2f713370ff34715264ef24a8176cb1caa53c28d Mon Sep 17 00:00:00 2001 From: Milan Date: Fri, 5 Apr 2024 00:35:15 -0400 Subject: [PATCH 13/14] julia> ] space for Navid! --- docs/src/installation.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/installation.md b/docs/src/installation.md index 341e71b52..ed0995744 100644 --- a/docs/src/installation.md +++ b/docs/src/installation.md @@ -8,7 +8,7 @@ julia> Pkg.add("SpeedyWeather") ``` or, equivalently, (`]` opens the package manager) ```julia -julia>] add SpeedyWeather +julia> ] add SpeedyWeather ``` which will automatically install the [latest release](https://github.com/SpeedyWeather/SpeedyWeather.jl/releases) and all necessary dependencies. If you run into any troubles please raise an @@ -20,7 +20,7 @@ julia> Pkg.add(url="https://github.com/SpeedyWeather/SpeedyWeather.jl", rev="mai ``` In a similar manner, other branches can be also installed. You can also type, equivalently, ```julia -julia>] add https://github.com/SpeedyWeather/SpeedyWeather.jl#main +julia> ] add https://github.com/SpeedyWeather/SpeedyWeather.jl#main ``` ## Compatibility with Julia versions From 75d3d0eaa3567495eeb24b0c841e9980fba70843 Mon Sep 17 00:00:00 2001 From: Milan Date: Fri, 5 Apr 2024 00:49:31 -0400 Subject: [PATCH 14/14] swap () and [] --- docs/src/speedytransforms.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/speedytransforms.md b/docs/src/speedytransforms.md index 730ed6662..c46a796b2 100644 --- a/docs/src/speedytransforms.md +++ b/docs/src/speedytransforms.md @@ -241,7 +241,7 @@ heatmap(map, title="k⁻²-distributed noise") save("random_noise.png", ans) # hide nothing # hide ``` -!(Random noise)[random_noise.png] +![Random noise](random_noise.png) You can always access the underlying data in `map` via `map.data` in case you need to get rid of the wrapping into a grid again! @@ -412,7 +412,7 @@ heatmap(η_grid, title="Geostrophic interface displacement η [m]") save("eta_geostrophic.png", ans) # hide nothing # hide ``` -!(Geostrophic eta)[eta_geostrophic.png] +![Geostrophic eta](eta_geostrophic.png) Which is the interface displacement assuming geostrophy. The actual interface displacement contains also ageostrophy @@ -422,7 +422,7 @@ heatmap(η_grid2, title="Interface displacement η [m] with ageostrophy") save("eta_ageostrophic.png", ans) # hide nothing # hide ``` -!(Ageostrophic eta)[eta_ageostrophic.png] +![Ageostrophic eta](eta_ageostrophic.png) Strikingly similar! The remaining differences are the ageostrophic motions but also note that the mean is off. This is because geostrophy only use/defines the gradient