From 35314c04e5930ad4eb41f4e2a6d4663227439c5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= Date: Fri, 5 Apr 2024 01:09:03 -0400 Subject: [PATCH] Held-Suarez example, NoTypes for surface fluxes (#514) * Split into 2D/3D examples * NoTypes for surface fluxes * Held-Suarez example added * SpeedyTransforms docs with Makie too --- README.md | 2 +- docs/make.jl | 3 +- docs/src/convection.md | 29 ++- docs/src/{examples.md => examples_2D.md} | 124 ++----------- docs/src/examples_3D.md | 219 +++++++++++++++++++++++ docs/src/forcing_drag.md | 2 +- docs/src/how_to_run_speedy.md | 4 +- docs/src/installation.md | 6 +- docs/src/large_scale_condensation.md | 26 ++- docs/src/output.md | 2 +- docs/src/radiation.md | 3 +- docs/src/speedytransforms.md | 44 +++-- src/models/primitive_dry.jl | 7 +- src/models/primitive_wet.jl | 12 +- src/physics/surface_fluxes.jl | 96 +++++----- src/physics/temperature_relaxation.jl | 16 +- 16 files changed, 405 insertions(+), 190 deletions(-) rename docs/src/{examples.md => examples_2D.md} (72%) create mode 100644 docs/src/examples_3D.md 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/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..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 +## 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 @@ -152,7 +162,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 @@ -182,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/examples.md b/docs/src/examples_2D.md similarity index 72% rename from docs/src/examples.md rename to docs/src/examples_2D.md index 4b90d3796..601b22293 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 @@ -137,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) @@ -152,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. @@ -160,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 @@ -176,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) @@ -304,7 +309,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 +317,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..1f0293c8c --- /dev/null +++ b/docs/src/examples_3D.md @@ -0,0 +1,219 @@ +# 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 + +```@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), + boundary_layer_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_heat_flux = NoSurfaceHeatFlux(), + + # 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 +(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. + +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 + +```@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) +model.feedback.verbose = false # hide +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) +model.feedback.verbose = false # hide +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. 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) 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..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 @@ -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. 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) 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`. 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) ``` diff --git a/docs/src/speedytransforms.md b/docs/src/speedytransforms.md index 71b839336..c46a796b2 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 diff --git a/src/models/primitive_dry.jl b/src/models/primitive_dry.jl index 6b2fb10b6..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<:AbstractSurfaceHeat, + 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) @@ -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..7a4c5eda8 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<:AbstractSurfaceHeatFlux, + EV<:AbstractSurfaceEvaporation, LSC<:AbstractCondensation, CV<:AbstractConvection, SW<:AbstractShortwave, @@ -82,8 +82,8 @@ 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) - evaporation::EV = SurfaceEvaporation(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) 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..d8f78c36e 100644 --- a/src/physics/surface_fluxes.jl +++ b/src/physics/surface_fluxes.jl @@ -1,28 +1,30 @@ -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 AbstractSurfaceHeatFlux <: 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) + surface_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() +initialize!(::SurfaceThermodynamicsConstant,::PrimitiveEquation) = nothing function surface_thermodynamics!( column::ColumnVariables, ::SurfaceThermodynamicsConstant, - atmosphere::AbstractAtmosphere, model::PrimitiveWet) # surface value is same as lowest model level @@ -30,21 +32,27 @@ 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 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) 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 -export SurfaceSensibleHeat -Base.@kwdef struct SurfaceSensibleHeat{NF<:AbstractFloat} <: AbstractSurfaceHeat +## SENSIBLE HEAT FLUX +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 @@ -122,14 +139,15 @@ Base.@kwdef struct SurfaceSensibleHeat{NF<:AbstractFloat} <: AbstractSurfaceHeat max_flux::NF = 100 end -SurfaceSensibleHeat(SG::SpectralGrid; kwargs...) = SurfaceSensibleHeat{SG.NF}(; kwargs...) +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, - atmosphere::AbstractAtmosphere + heat_flux::SurfaceHeatFlux, + model::PrimitiveEquation, ) - 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 <: AbstractSurfaceHeatFlux 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 @@ -191,22 +216,11 @@ Base.@kwdef struct SurfaceEvaporation{NF<:AbstractFloat} <: AbstractEvaporation end SurfaceEvaporation(SG::SpectralGrid; kwargs...) = SurfaceEvaporation{SG.NF}(; kwargs...) +initialize!(::SurfaceEvaporation, ::PrimitiveWet) = nothing -# 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 +228,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 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