Analysing a simulation
While you can analyze a SpeedyWeather simulation through its NetCDF output (i.e., offline analysis), as most users will be used to with other models, you can also reuse a lot of functionality from SpeedyWeather interactively for analysis. This makes SpeedyWeather beyond being a model for the atmospheric general circulation also a library with many functions for the analysis of simulations. Often this also avoids the two-language problem that you will face if you run a simulation with a model in one language but then do the data analysis in another, treating the model as a blackbox although it likely has many of the functions you will need for analysis already defined. With SpeedyWeather we try to avoid this and are working towards a more unified approach in atmospheric modelling where simulation and analysis are done interactively with the same library: SpeedyWeather.jl.
Advantages of online analysis
Now you could run a SpeedyWeather simulation, and analyse the NetCDF output but that comes with several issues related to accuracy
- If you use a reduced grid for the simulation, then the output will (by default) be interpolated on a full grid. This interpolation introduces an error.
- Computing integrals over gridded data on the sphere by weighting every grid point according to its area is not the most accurate numerical integration.
- Computing gradients over gridded data comes with similar issues. While our RingGrids are always equidistant in longitude, they are not necessarily in latitude.
The first point you can avoid by running a simulation on one of the full grids that are implemented, see SpectralGrid. But that also impacts the simulation and for various reasons we don't run on full grids by default.
The second point you can address by defining a more advanced numerical integration scheme, but that likely requires you to depend on external libraries and then, well, you could also just depend on SpeedyWeather.jl directly, because we have to do these computations internally anyway. Similar for the third point, many gradients have to be computed on every time step and we do that with spectral transforms to reduce the discretization error.
The following contains a (hopefully growing) list of examples of how a simulation can be analysed interactively. We call this online analysis because you are directly using the functionality from the model as if it was a library. For simplicity, we use the shallow water model to demonstrate this.
Mass conservation
In the absence of sources and sinks for the interface displacement $\eta$ in the shallow water equations, total mass (or equivalently volume as density is constant) is conserved. The total volume is defined as the integral of the dynamic layer thickness $h = \eta + H - H_b$ ($H$ is the layer thickness at rest, $H_b$ is orography) over the surface $A$ of the sphere
\[M = \iint h dA = \iint \left(\eta + H - H_b\right) dA\]
to check for conservation we want to assess that
\[\frac{\partial M}{\partial t} = 0\]
And because $V = \iint H dA - \iint H_b dA$, the total volume at rest, is a constant ($H$ is a global constant, the orography $H_b$ does not change with time) we can just check whether $\iint \eta dA$ changes over time. Instead of computing this integral in grid-point space, we use the spectral $\eta$ whereby the coefficient of the first spherical harmonic (the $l = m = 0$ mode, or wavenumber 0, see Spherical Harmonic Transform) encodes the global average.
using SpeedyWeather
+spectral_grid = SpectralGrid(trunc=31, nlayers=1)
+model = ShallowWaterModel(;spectral_grid)
+simulation = initialize!(model)Simulation{ShallowWaterModel}
+├ prognostic_variables::PrognosticVariables{...}
+├ diagnostic_variables::DiagnosticVariables{...}
+└ model::ShallowWaterModel{...}Now we check $\eta_{0,0}$ the $l = m = 0$ coefficent of the inital conditions of that simulation with
simulation.prognostic_variables.pres[1][1]1.1815667f0 + 0.0f0im[1][1] pulls the first Leapfrog time step and of that the first element of the underlying LowerTriangularMatrix which is the coefficient of the $l = m = 0$ harmonic. Its imaginary part is always zero (which is true for any zonal harmonic $m=0$ as its imaginary part would just unnecessarily rotate something zonally constant in zonal direction), so you can real it. Also for spherical harmonic transforms there is a norm of the sphere by which you have to divide to get your mean value in the original units
a = model.spectral_transform.norm_sphere # = 2√π = 3.5449078
+η_mean = real(simulation.prognostic_variables.pres[1][1]) / a0.3333138f0So the initial conditions in this simulation are such that the global mean interface displacement is that value in meters. You would need to multiply by the area of the sphere $4\pi R^2$ (radius $R$) to get the actual integral from above, but because that doesn't change with time either, we just want to check that η_mean doesn't change with time. Which is equivalent to $\partial_t \iint \eta dA = 0$ and so volume conservation and because density is constant also mass conservation. Let's check what happens after the simulation ran for some days
run!(simulation, period=Day(10))
+
+# now we check η_mean again
+η_mean_later = real(simulation.prognostic_variables.pres[1][1]) / a0.3333138f0which is exactly the same. So mass is conserved, woohoo.
Insight from a numerical perspective: The tendency of $\eta$ is $\partial_t \eta = -\nabla \cdot (\mathbf{u} h)$ which is a divergence of a flux. Calculating the divergence in spherical harmonics always sets the $l=m=0$ mode to zero exactly (the gradient of a periodic variable has zero mean) so the time integration here is always with an exactly zero tendency.
Energy
The total energy in the shallow water equations is the sum of the kinetic energy density $\frac{1}{2}(u^2 + v^2)$ and the potential energy density $gz$ integrated over the total volume (times $h=\eta+H-H_b$ for the vertical then integrated over the sphere $\iint dA$).
\[\begin{align} +E & = \iint \left[ \int_{H_b}^{\eta}\frac{1}{2}\left(u^2 + v^2 + gz\right)dz \right]dA \\ + & = \iint \frac{1}{2}\left(u^2 + v^2 + gh\right)h dA +\end{align}\]
In contrast to the Mass conservation which, with respect to the spectral transform, is a linear calculation, here we need to multiply variables, which has to be done in grid-point space. Then we can transform to spectral space for the global integral as before. Let us define a total_energy function as
using SpeedyWeather
+function total_energy(u, v, η, model)
+ H = model.atmosphere.layer_thickness
+ Hb = model.orography.orography
+ g = model.planet.gravity
+
+ h = @. η + H - Hb # layer thickness between the bottom and free surface
+ E = @. h/2*(u^2 + v^2) + g*h^2 # vertically-integrated mechanical energy
+
+ # transform to spectral, take l=m=0 mode at [1] and normalize for mean
+ return E_mean = real(transform(E)[1]) / model.spectral_transform.norm_sphere
+endtotal_energy (generic function with 1 method)So at the current state of our simulation we have a total energy (per square meter as we haven't multiplied by the surface area of the planet).
# flat copies for convenience
+u = simulation.diagnostic_variables.grid.u_grid[:, 1]
+v = simulation.diagnostic_variables.grid.v_grid[:, 1]
+η = simulation.diagnostic_variables.grid.pres_grid
+
+TE = total_energy(u, v, η, model)6.750372f8with units of $m^3 s^{-2}$ (multiplying by surface area of the sphere and density of the fluid would turn it into joule = $kg \, m^2 s^{-2}$). To know in general where to find the respective variables $u, v, \eta$ inside our simulation object see Prognostic variables and Diagnostic variables. Now let us continue the simulation
run!(simulation, period=Day(10))
+
+# we don't need to reassign u, v, η as they were flat copies
+# pointing directly to the diagnostic variables inside the simulation
+# which got updated during run!
+TE_later = total_energy(u, v, η, model)6.746546f8So the total energy has somewhat changed, it decreased to
TE_later/TE0.9994332f0of its previous value over 10 days.
While technically energy should be conserved in an unforced system, numerically this is rarely exactly the case. We need some Horizontal diffusion for numerical stability and also the time integration is dissipative due to temporal filtering, see Time integration. Note that the energy here is inversely cascading to larger scales, and this makes the dissipation of energy through Horizontal diffusion really inefficient, because in spectral space, standard Laplacian diffusion is proportional to $k^2$ where $k$ is the wavenumber. By default, we use a 4th power Laplacian, so the diffusion here is proportional to $k^8$.
Potential vorticity
Potential vorticity in the shallow water equations is defined as
\[q = \frac{f + \zeta}{h}\]
with $f$ the Coriolis parameter, $\zeta$ the relative vorticity, and $h$ the layer thickness as before. We can calculate this conveniently directly on the model grid (whichever you chose) as
# vorticity
+ζ = simulation.diagnostic_variables.grid.vor_grid[:,1]
+f = coriolis(ζ) # create f on that grid
+
+# layer thickness
+η = simulation.diagnostic_variables.grid.pres_grid
+H = model.atmosphere.layer_thickness
+Hb = model.orography.orography
+h = @. η + H - Hb
+
+# potential vorticity
+q = @. (f + ζ) / hand we can compare the relative vorticity field to
using CairoMakie
+heatmap(ζ, title="Relative vorticity [1/s]")
the potential vorticity
heatmap(q, title="Potential vorticity [1/m/s]")
Absolute angular momentum
Similar to the total mass, in the absence of sources and sinks for momentum, total absolute angular momentum (AAM) defined as
\[\Lambda = \iint \left(ur + \Omega r^2\right)h dA\]
should be conserved ($\partial_t\Lambda = 0$). Here $u$ is the zonal velocity, $\Omega$ the angular velocity of the Earth, $r = R \cos\phi$ the momentum arm at latitude $\phi$, and $R$ the radius of Earth.
Following previous examples, let us define a total_angular_momentum function as
using SpeedyWeather
+
+function total_angular_momentum(u, η, model)
+ H = model.atmosphere.layer_thickness
+ Hb = model.orography.orography
+ R = model.spectral_grid.radius
+ Ω = model.planet.rotation
+
+ r = R * cos.(model.geometry.lats) # momentum arm for every grid point
+
+ h = @. η + H - Hb # layer thickness between the bottom and free surface
+ Λ = @. (u*r + Ω*r^2) * h # vertically-integrated AAM
+
+ # transform to spectral, take l=m=0 mode at [1] and normalize for mean
+ return Λ_mean = real(transform(Λ)[1]) / model.spectral_transform.norm_sphere
+endtotal_angular_momentum (generic function with 1 method)Anytime we stop the simulation, we can calculate $\Lambda$ using this function (ignoring the multiplication by $4\pi R^2$ to get total $\Lambda$).
# use u, η from current state of simulation
+Λ_current = total_angular_momentum(u, η, model)1.65041281411102e13So after some days of integration, we would get another $\Lambda$ with
run!(simulation, period=Day(10))
+
+# u, η got updated during run!
+Λ_later = total_angular_momentum(u, η, model)
+Λ_later / Λ_current0.9981800468569637and measure its relative change. Similar to energy, in the numerical integration, Λ is not exactly conserved due to Horizontal diffusion.
Circulation
Total circulation is defined as the area-integrated absolute vorticity:
\[C = \iint \left(\zeta + f\right) dA\]
Following previous fashion, we define a function total_circulation for this
function total_circulation(ζ, model)
+ f = coriolis(ζ) # create f on the grid of ζ
+ C = ζ .+ f # absolute vorticity
+ # transform to spectral, take l=m=0 mode at [1] and normalize for mean
+ return C_mean = real(transform(C)[1]) / model.spectral_transform.norm_sphere
+end
+
+total_circulation(ζ, model)4.321589f-13Note that the area integral of relative vorticity $\zeta$ and planetary vorticity $f$ over the whole surface of a sphere are analytically exactly zero. Numerically, C_mean should be a small number but may not be exactly zero due to numerical precision and errors in the spectral transform.
Potential enstrophy
The total potential enstrophy is defined as the second-moment of potential vorticity $q$
\[Q = \iint \frac{1}{2}q^2 dA\]
In the absence of source and sink for potential vorticiy, this quantity should also conserve during the integration.
We define a function $total_enstrophy$ for this
function total_enstrophy(ζ, η, model)
+ # constants from model
+ H = model.atmosphere.layer_thickness
+ Hb = model.orography.orography
+ f = coriolis(ζ) # create f on the grid
+
+ h = @. η + H - Hb # thickness
+ q = @. (ζ + f) / h # Potential vorticity
+ Q = @. q^2 / 2 # Potential enstrophy
+
+ # transform to spectral, take l=m=0 mode at [1] and normalize for mean
+ return Q_mean = real(transform(Q)[1]) / model.spectral_transform.norm_sphere
+endtotal_enstrophy (generic function with 1 method)Then by evaluating Q_mean at different time steps, one can similarly check how $Q$ is changing over time.
Q = total_enstrophy(ζ, η, model)
+run!(simulation, period=Day(10))
+Q_later = total_enstrophy(ζ, η, model)
+Q_later/Q0.98658925f0Note that the turbulent nature of the shallow water model (or generally 2D turbulence) cascades enstrophy to smaller scales where it is removed by Horizontal diffusion for numerical stability. As a result, it is decreasing more quickly than energy.
Online diagnostics
Now we want to calculate all the above global diagnostics periodically during a simulation. For that we will use Callbacks, which let us inject code into a simulation that is executed after every time step (or at any other scheduled time, see Schedules).
So we define a function global_diagnostics to calculate the integrals together. We could reuse the functions like total_enstrophy from above but we also want to show how to global integral $\iint dV$ can be written more efficiently
# define a global integral, reusing a precomputed SpectralTransform S
+# times surface area of sphere omitted
+function ∬dA(v, h, S::SpectralTransform)
+ return real(transform(v .* h, S)[1]) / S.norm_sphere
+end
+
+# use SpectralTransform from model
+∬dA(v, h, model::AbstractModel) = ∬dA(v, h, model.spectral_transform)∬dA (generic function with 2 methods)By reusing model.spectral_transform we do not have to re-precompute the spectral tranform on every call to transform. Providing the spectral transform from model as the 2nd argument simply reuses a previously precomputed spectral transform which is much faster and uses less memory.
Now the global_diagnostics function is defined as
function global_diagnostics(u, v, ζ, η, model)
+
+ # constants from model
+ NF = model.spectral_grid.NF # number format used
+ H = model.atmosphere.layer_thickness
+ Hb = model.orography.orography
+ R = model.spectral_grid.radius
+ Ω = model.planet.rotation
+ g = model.planet.gravity
+
+ r = NF.(R * cos.(model.geometry.lats)) # create r on that grid
+ f = coriolis(u) # create f on that grid
+
+ h = @. η + H - Hb # thickness
+ q = @. (ζ + f) / h # potential vorticity
+ λ = @. u * r + Ω * r^2 # angular momentum (in the right number format NF)
+ k = @. (u^2 + v^2) / 2 # kinetic energy
+ p = @. g * h / 2 # potential energy
+ z = @. q^2 / 2 # potential enstrophy
+
+ M = ∬dA(1, h, model) # mean mass
+ C = ∬dA(q, h, model) # mean circulation
+ Λ = ∬dA(λ, h, model) # mean angular momentum
+ K = ∬dA(k, h, model) # mean kinetic energy
+ P = ∬dA(p, h, model) # mean potential energy
+ Q = ∬dA(z, h, model) # mean potential enstrophy
+
+ return M, C, Λ, K, P, Q
+end
+
+# unpack diagnostic variables and call global_diagnostics from above
+function global_diagnostics(diagn::DiagnosticVariables, model::AbstractModel)
+ u = diagn.grid.u_grid[:, 1]
+ v = diagn.grid.v_grid[:, 1]
+ ζR = diagn.grid.vor_grid[:, 1]
+ η = diagn.grid.pres_grid
+
+ # vorticity during simulation is scaled by radius R, unscale here
+ ζ = ζR ./ diagn.scale[]
+
+ return global_diagnostics(u, v, ζ, η, model)
+endglobal_diagnostics (generic function with 2 methods)The prognostic variables vorticity and divergence are scaled with the radius of the sphere during a simulation, see Radius scaling. This did not apply above because we only analyzed vorticity before or after the simulation, i.e. outside of the run!(simulation) call. The radius scaling is only applied just before the time integration and is undone directly after it. However, because now we are accessing the vorticity during the simulation we need to unscale the vorticity (and divergence) manually. General recommendation is to divide by diagn.scale[] (and not radius) as diagn.scale[] always reflects whether a vorticity and divergence are currently scaled (scale = radius) or not (scale = 1).
Then we define a new callback GlobalDiagnostics subtype of SpeedyWeather's AbstractCallback and define new methods of initialize!, callback! and finalize! for it (see Callbacks for more details)
# define a GlobalDiagnostics callback and the fields it needs
+Base.@kwdef mutable struct GlobalDiagnostics <: SpeedyWeather.AbstractCallback
+ timestep_counter::Int = 0
+
+ time::Vector{DateTime} = []
+ M::Vector{Float64} = [] # mean mass per time step
+ C::Vector{Float64} = [] # mean circulation per time step
+ Λ::Vector{Float64} = [] # mean angular momentum per time step
+ K::Vector{Float64} = [] # mean kinetic energy per time step
+ P::Vector{Float64} = [] # mean potential energy per time step
+ Q::Vector{Float64} = [] # mean enstrophy per time step
+end
+
+# define how to initialize a GlobalDiagnostics callback
+function SpeedyWeather.initialize!(
+ callback::GlobalDiagnostics,
+ progn::PrognosticVariables,
+ diagn::DiagnosticVariables,
+ model::AbstractModel,
+)
+ # replace with vector of correct length
+ n = progn.clock.n_timesteps + 1 # +1 for initial conditions
+ callback.time = zeros(DateTime, n)
+ callback.M = zeros(n)
+ callback.C = zeros(n)
+ callback.Λ = zeros(n)
+ callback.K = zeros(n)
+ callback.P = zeros(n)
+ callback.Q = zeros(n)
+
+ M, C, Λ, K, P, Q = global_diagnostics(diagn, model)
+
+ callback.time[1] = progn.clock.time
+ callback.M[1] = M # set initial conditions
+ callback.C[1] = C # set initial conditions
+ callback.Λ[1] = Λ # set initial conditions
+ callback.K[1] = K # set initial conditions
+ callback.P[1] = P # set initial conditions
+ callback.Q[1] = Q # set initial conditions
+
+ callback.timestep_counter = 1 # (re)set counter to 1
+
+ return nothing
+end
+
+# define what a GlobalDiagnostics callback does on every time step
+function SpeedyWeather.callback!(
+ callback::GlobalDiagnostics,
+ progn::PrognosticVariables,
+ diagn::DiagnosticVariables,
+ model::AbstractModel,
+)
+ callback.timestep_counter += 1
+ i = callback.timestep_counter
+
+ M, C, Λ, K, P, Q = global_diagnostics(diagn, model)
+
+ # store current time and diagnostics for timestep i
+ callback.time[i] = progn.clock.time
+ callback.M[i] = M
+ callback.C[i] = C
+ callback.Λ[i] = Λ
+ callback.K[i] = K
+ callback.P[i] = P
+ callback.Q[i] = Q
+end
+
+using NCDatasets
+
+# define how to finalize a GlobalDiagnostics callback after simulation finished
+function SpeedyWeather.finalize!(
+ callback::GlobalDiagnostics,
+ progn::PrognosticVariables,
+ diagn::DiagnosticVariables,
+ model::AbstractModel,
+)
+ n_timesteps = callback.timestep_counter
+
+ # create a netCDF file in current path
+ ds = NCDataset(joinpath(pwd(), "global_diagnostics.nc"), "c")
+
+ # save diagnostics variables within
+ defDim(ds, "time", n_timesteps)
+ defVar(ds, "time", callback.time, ("time",))
+ defVar(ds, "mass", callback.M, ("time",))
+ defVar(ds, "circulation", callback.C, ("time",))
+ defVar(ds, "angular momentum", callback.Λ, ("time",))
+ defVar(ds, "kinetic energy", callback.K, ("time",))
+ defVar(ds, "potential energy", callback.P, ("time",))
+ defVar(ds, "potential enstrophy", callback.Q, ("time",))
+
+ close(ds)
+
+ return nothing
+endNote that callback! will execute every time step. If execution is only desired periodically, you can use Schedules. At finalize! we decide to write the timeseries of our global diagnostics as netCDF file via NCDatasets.jl to the current path pwd(). We need to add using NCDatasets here, as SpeedyWeather does not re-export the functionality therein.
Now we create a GlobalDiagnostics callback, add it to the model with key :global_diagnostics (you get a random key if not provided) and reinitialize the simulation to start from the initial conditions.
# don't name it global_diagnostics because that's a function already!
+diagnostics_recorder = GlobalDiagnostics()
+add!(model.callbacks, :diagnostics_recorder => diagnostics_recorder)
+simulation = initialize!(model)
+
+run!(simulation, period=Day(20))Then one could check the output file global_diagnostics.nc, or directly use the callback through its key :diagnostics_recorder as we do here
using CairoMakie
+
+# unpack callback
+(; M, C, Λ, K, P, Q) = model.callbacks[:diagnostics_recorder]
+t = model.callbacks[:diagnostics_recorder].time
+days = [Second(ti - t[1]).value/3600/24 for ti in t]
+
+fig = Figure();
+axs = [Axis(fig[row, col]) for row in 1:3, col in 1:2]
+
+lines!(axs[1,1], days, M)
+axs[1,1].title = "Mass"
+hidexdecorations!(axs[1,1])
+
+lines!(axs[2,1], days, Λ)
+axs[2,1].title = "Angular momentum"
+hidexdecorations!(axs[2,1])
+
+lines!(axs[3,1], days, P)
+axs[3,1].title = "Potential energy"
+axs[3,1].xlabel = "time [day]"
+
+lines!(axs[1,2], days, C)
+axs[1,2].title = "Circulation"
+hidexdecorations!(axs[1,2])
+
+lines!(axs[2,2], days, K)
+axs[2,2].title = "Kinetic energy"
+hidexdecorations!(axs[2,2])
+
+lines!(axs[3,2], days, Q)
+axs[3,2].title = "Potential enstrophy"
+axs[3,2].xlabel = "time [day]"
+fig
![](\"https://upload.wikimedia.org/wikipedia/commons/1/12/Rotating_spherical_harmonics.gif\")
","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Every row represents an order l geq 0, starting from l=0 at the top. Every column represents an order m geq 0, starting from m=0 on the left. The coefficients of these spherical harmonics are directly mapped into a matrix a_lm as ","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":" m \nl a_00 \n a_10 a_11 \n a_20 a_12 a_22 \n a_30 a_13 a_23 a_33","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"which is consistently extended for higher degrees and orders. Consequently, all spectral fields are lower-triangular matrices with complex entries. The upper triangle excluding the diagonal are zero. Note that internally vector fields include an additional degree, such that l_max = m_max + 1 (see Derivatives in spherical coordinates for more information). The harmonics with a_l0 (the first column) are also called zonal harmonics as they are constant with longitude phi. The harmonics with a_ll (the main diagonal) are also called sectoral harmonics as they essentially split the sphere into 2l sectors in longitude phi without a zero-crossing in latitude.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"For correctness it is mentioned here that SpeedyWeather.jl uses a LowerTriangularMatrix type to store the spherical harmonic coefficients. By doing so, the upper triangle is actually not explicitly stored and the data technically unravelled into a vector, but this is hidden as much as possible from the user. For more details see LowerTriangularMatrices.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"info: Array indices\nFor a spectral field a note that due to Julia's 1-based indexing the coefficient a_lm is obtained via a[l+1, m+1]. Alternatively, we may index over 1-based l, m but a comment is usually added for clarification.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Fortran SPEEDY does not use the same spectral packing as SpeedyWeather.jl. The alternative packing l m therein uses l=m and m=l-m as summarized in the following table.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"degree l order m l=m m=l-m\n0 0 0 0\n1 0 0 1\n1 1 1 0\n2 0 0 2\n2 1 1 1\n2 2 2 0\n3 0 0 3\n... ... ... ...","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"This alternative packing uses the top-left triangle of a coefficient matrix, and the degrees and orders from above are stored at the following indices","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":" m \nl a_00 a_10 a_20 a_30\n a_11 a_21 a_31 \n a_22 a_32 \n a_33 ","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"This spectral packing is not used in SpeedyWeather.jl but illustrated here for completeness and comparison with Fortran SPEEDY.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"SpeedyWeather.jl uses triangular truncation such that only spherical harmonics with l leq l_max and m leq m_max are explicitly represented. This is usually described as Tm_max, with l_max = m_max (although in vector quantities require one more degree l in the recursion relation of meridional gradients). For example, T31 is the spectral resolution with l_max = m_max = 31. Note that the degree l and order m are mathematically 0-based, such that the corresponding coefficient matrix is of size 32x32.","category":"page"},{"location":"spectral_transform/#Available-horizontal-resolutions","page":"Spherical Harmonic Transform","title":"Available horizontal resolutions","text":"","category":"section"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Technically, SpeedyWeather.jl supports arbitrarily chosen resolution parameter trunc when creating the SpectralGrid that refers to the highest non-zero degree l_max that is resolved in spectral space. SpeedyWeather.jl will always try to choose an easily-Fourier transformable[FFT] size of the grid, but as we use FFTW.jl there is quite some flexibility without performance sacrifice. However, this has traditionally lead to typical resolutions that we also use for testing we therefore recommend to use. They are as follows with more details below","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"trunc nlon nlat Delta x\n31 (default) 96 48 400 km\n42 128 64 312 km\n63 192 96 216 km\n85 256 128 165 km\n127 384 192 112 km\n170 512 256 85 km\n255 768 384 58 km\n341 1024 512 43 km\n511 1536 768 29 km\n682 2048 1024 22 km\n1024 3072 1536 14 km\n1365 4092 2048 11 km","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Some remarks on this table","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"This assumes the default quadratic truncation, you can always adapt the grid resolution via the dealiasing option, see Matching spectral and grid resolution\nnlat refers to the total number of latitude rings, see Grids. With non-Gaussian grids, nlat will be one one less, e.g. 47 instead of 48 rings.\nnlon is the number of longitude points on the Full Gaussian Grid, for other grids there will be at most these number of points around the Equator.\nDelta x is the horizontal resolution. For a spectral model there are many ways of estimating this[Randall2021]. We use here the square root of the average area a grid cell covers, see Effective grid resolution","category":"page"},{"location":"spectral_transform/#Effective-grid-resolution","page":"Spherical Harmonic Transform","title":"Effective grid resolution","text":"","category":"section"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"There are many ways to estimate the effective grid resolution of spectral models[Randall2021]. Some of them are based on the wavelength a given spectral resolution allows to represent, others on the total number of real variables per area. However, as many atmospheric models do represent a considerable amount of physics on the grid (see Parameterizations) there is also a good argument to include the actual grid resolution into this estimate and not just the spectral resolution. We therefore use the average grid cell area to estimate the resolution","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Delta x = sqrtfrac4pi R^2N","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"with N number of grid points over a sphere with radius R. However, we have to acknowledge that this usually gives higher resolution compared to other methods of estimating the effective resolution, see [Randall2021] for a discussion. You may therefore need to be careful to make claims that, e.g. trunc=85 can resolve the atmospheric dynamics at a scale of 165km.","category":"page"},{"location":"spectral_transform/#Derivatives-in-spherical-coordinates","page":"Spherical Harmonic Transform","title":"Derivatives in spherical coordinates","text":"","category":"section"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Horizontal gradients in spherical coordinates are defined for a scalar field A and the latitudes theta and longitudes lambda as","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"nabla A = left(frac1Rcosthetafracpartial Apartial lambda frac1Rfracpartial Apartial theta right)","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"However, the divergence of a vector field mathbfu = (u v) includes additional cos(theta) scalings","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"nabla cdot mathbfu = frac1Rcosthetafracpartial upartial lambda +\nfrac1Rcosthetafracpartial (v costheta)partial theta","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"and similar for the curl","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"nabla times mathbfu = frac1Rcosthetafracpartial vpartial lambda -\nfrac1Rcosthetafracpartial (u costheta)partial theta","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"The radius of the sphere (i.e. Earth) is R. The zonal gradient scales with 1cos(theta) as the longitudes converge towards the poles (note that theta describes latitudes here, definitions using colatitudes replace the cos with a sin.)","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Starting with a spectral field of vorticity zeta and divergence mathcalD one can obtain stream function Psi and velocity potential Phi by inverting the Laplace operator nabla^2:","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Psi = nabla^-2zeta quad Phi = nabla^-2mathcalD","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"The velocities u v are then obtained from (u v) = nabla^botPsi + nablaPhi following the definition from above and nabla^bot = (-R^-1partial_theta (Rcostheta)^-1partial_lambda)","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"beginaligned\nu = -frac1Rpartial_thetaPsi + frac1Rcosthetapartial_lambdaPhi \nv = +frac1Rpartial_thetaPhi + frac1Rcosthetapartial_lambdaPsi\nendaligned","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"How the operators nabla nabla times nabla cdot can be implemented with spherical harmonics is presented in the following sections. However, note that the actually implemented operators differ slightly in their scaling with respect to the radius R and the cosine of latitude cos(theta). For further details see Gradient operators which describes those as implemented in the SpeedyTransforms module. Also note that the equations in SpeedyWeather.jl are scaled with the radius R^2 (see Radius scaling) which turns most operators into non-dimensional operators on the unit sphere anyway.","category":"page"},{"location":"spectral_transform/#Zonal-derivative","page":"Spherical Harmonic Transform","title":"Zonal derivative","text":"","category":"section"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"The zonal derivative of a scalar field Psi in spectral space is the zonal derivative of all its respective spherical harmonics Psi_lm(phi theta) (now we use phi for longitudes to avoid confusion with the Legendre polynomials lambda_lm)","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"v_lm = frac1R cos(theta) fracpartialpartial phi left( lambda_l^m(costheta) e^imphi right) =\nfracimR cos(theta) lambda_l^m(costheta) e^imphi = fracimR cos(theta) Psi_lm","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"So for every spectral harmonic, cos(theta)v_lm is obtained from Psi_lm via a multiplication with imR. Unscaling the cos(theta)-factor is done after transforming the spectral coefficients v_lm into grid-point space. As discussed in Radius scaling, SpeedyWeather.jl scales the stream function as tildePsi = R^-1Psi such that the division by radius R in the gradients can be omitted. The zonal derivative becomes therefore effectively for each spherical harmonic a scaling with its (imaginary) order im. The spherical harmonics are essentially just a Fourier transform in zonal direction and the derivative a multiplication with the respective wave number m times imaginary i.","category":"page"},{"location":"spectral_transform/#Meridional-derivative","page":"Spherical Harmonic Transform","title":"Meridional derivative","text":"","category":"section"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"The meridional derivative of the spherical harmonics is a derivative of the Legendre polynomials for which the following recursion relation applies[Randall2021], [Durran2010], [GFDL], [Orszag70]","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"costheta fracdP_l mdtheta = -lepsilon_l+1 mP_l+1 m + (l+1)epsilon_l mP_l-1 m","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"with recursion factors","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"epsilon_l m = sqrtfracl^2-m^24l^2-1","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"In the following we use the example of obtaining the zonal velocity u from the stream function Psi, which is through the negative meridional gradient. For the meridional derivative itself the leading minus sign has to be omitted. Starting with the spectral expansion","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Psi(lambda theta) = sum_l mPsi_l mP_l m(sintheta)e^imlambda","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"we multiply with -R^-1costhetapartial_theta to obtain","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"costhetaleft(-frac1Rpartial_thetaPsi right) = -frac1Rsum_l mPsi_l me^imlambdacosthetapartial_theta P_l m","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"at which point the recursion from above can be applied. Collecting terms proportional to P_l m then yields","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"(cos(theta)u)_l m = -frac1R(-(l-1)epsilon_l mPsi_l-1 m + (l+2)epsilon_l+1 mPsi_l+1 m)","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"To obtain the coefficient of each spherical harmonic l m of the meridional gradient of a spectral field, two coefficients at l-1 m and l+1 m have to be combined. This means that the coefficient of a gradient ((costheta) u)_lm is a linear combination of the coefficients of one higher and one lower degree Psi_l+1 m Psi_l-1 m. As the coefficient Psi_lm with ml are zero, the sectoral harmonics (l=m) of the gradients are obtained from the first off-diagonal only. However, the l=l_max harmonics of the gradients require the l_max-1 as well as the l_max+1 harmonics. As a consequence vector quantities like velocity components u v require one more degree l than scalar quantities like vorticity[Bourke72]. However, for easier compatibility all spectral fields in SpeedyWeather.jl use one more degree l, but scalar quantities should not make use of it. Equivalently, the last degree l is set to zero before the time integration, which only advances scalar quantities.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"In SpeedyWeather.jl, vector quantities like u v use therefore one more meridional mode than scalar quantities such as vorticity zeta or stream function Psi. The meridional derivative in SpeedyWeather.jl also omits the 1R-scaling as explained for the Zonal derivative and in Radius scaling.","category":"page"},{"location":"spectral_transform/#Divergence-and-curl-in-spherical-harmonics","page":"Spherical Harmonic Transform","title":"Divergence and curl in spherical harmonics","text":"","category":"section"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"The meridional gradient as described above can be applied to scalars, such as Psi and Phi in the conversion to velocities (u v) = nabla^botPsi + nablaPhi, however, the operators curl nabla times and divergence nabla cdot in spherical coordinates involve a costheta scaling before the meridional gradient is applied. How to translate this to spectral coefficients has to be derived separately[Randall2021], [Durran2010].","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"The spectral transform of vorticity zeta is","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"zeta_l m = frac12piint_-tfracpi2^tfracpi2int_0^2pi zeta(lambda theta)\nP_l m(sintheta) e^imlambda dlambda costheta dtheta","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Given that Rzeta = cos^-1partial_lambda v - cos^-1partial_theta (u costheta), we therefore have to evaluate a meridional integral of the form","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"int P_l m frac1cos theta partial_theta(u costheta) cos theta dtheta","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"which can be solved through integration by parts. As ucostheta = 0 at theta = pm tfracpi2 only the integral","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"= -int partial_theta P_l m (u costheta) dtheta = -int costheta partial_theta P_l m\n(fracucostheta) costheta dtheta","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"remains. Inserting the recurrence relation from the Meridional derivative turns this into","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"= -int left(-l epsilon_l+1 mP_l+1 m + (l+1)epsilon_l m P_l-1 m right) (fracucostheta)\ncos theta dtheta","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"Now we expand (tfracucostheta) but only the l m harmonic will project ontoP_l m. Let u^* = ucos^-1theta v^* = vcos^-1theta we then have in total","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"beginaligned\nRzeta_l m = imv^*_l m + (l+1)epsilon_l mu^*_l-1 m - lepsilon_l+1 mu^*_l+1 m \nRD_l m = imu^*_l m - (l+1)epsilon_l mv^*_l-1 m + lepsilon_l+1 mv^*_l+1 m \nendaligned","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"And the divergence D is similar, but (u v) to (-v u). We have moved the scaling with the radius R directly into zeta D as further described in Radius scaling.","category":"page"},{"location":"spectral_transform/#Laplacian","page":"Spherical Harmonic Transform","title":"Laplacian","text":"","category":"section"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"The spectral Laplacian is easily applied to the coefficients Psi_lm of a spectral field as the spherical harmonics are eigenfunctions of the Laplace operator nabla^2 in spherical coordinates with eigenvalues -l(l+1) divided by the radius squared R^2, i.e. nabla^2 Psi becomes tfrac-l(l+1)R^2Psi_lm in spectral space. For example, vorticity zeta and stream function Psi are related by zeta = nabla^2Psi in the barotropic vorticity model. Hence, in spectral space this is equivalent for every spectral mode of degree l and order m to","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"zeta_l m = frac-l(l+1)R^2Psi_l m","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"This can be easily inverted to obtain the stream function Psi from vorticity zeta instead. In order to avoid division by zero, we set Psi_0 0 here, given that the stream function is only defined up to a constant anyway.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"beginaligned\nPsi_l m = fracR^2-l(l+1)zeta_l m quad foralll m 0 \nPsi_0 0 = 0\nendaligned","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"See also Horizontal diffusion and Normalization of diffusion.","category":"page"},{"location":"spectral_transform/#U,-V-from-vorticity-and-divergence","page":"Spherical Harmonic Transform","title":"U, V from vorticity and divergence","text":"","category":"section"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"After having discussed the zonal and meridional derivatives with spherical harmonics as well as the Laplace operator, we can derive the conversion from vorticity zeta and divergence D (which are prognostic variables) to U=ucostheta V=vcostheta. Both are linear operations that act either solely on a given harmonic (the zonal gradient and the Laplace operator) or are linear combinations between one lower and one higher degree l (the meridional gradient). It is therefore computationally more efficient to compute U V directly from zeta D instead of calculating stream function and velocity potential first. In total we have","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"beginaligned\nU_l m = -fraciml(l+1)(RD)_l m + fracepsilon_l+1 ml+1(Rzeta)_l+1 m -\nfracepsilon_l ml(Rzeta)_l-1 m \nV_l m = -fraciml(l+1)(Rzeta)_l m - fracepsilon_l+1 ml+1(RD)_l+1 m +\nfracepsilon_l ml(RD)_l-1 m \nendaligned","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"We have moved the scaling with the radius R directly into zeta D as further described in Radius scaling.","category":"page"},{"location":"spectral_transform/#References","page":"Spherical Harmonic Transform","title":"References","text":"","category":"section"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[Malardel2016]: Malardel S, Wedi N, Deconinck W, Diamantakis M, Kühnlein C, Mozdzynski G, Hamrud M, Smolarkiewicz P. A new grid for the IFS. ECMWF newsletter. 2016; 146(23-28):321. doi: 10.21957/zwdu9u5i","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[Gorski2004]: Górski, Hivon, Banday, Wandelt, Hansen, Reinecke, Bartelmann, 2004. HEALPix: A FRAMEWORK FOR HIGH-RESOLUTION DISCRETIZATION AND FAST ANALYSIS OF DATA DISTRIBUTED ON THE SPHERE, The Astrophysical Journal. doi:10.1086/427976","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[Willmert2020]: Justin Willmert, 2020. justinwillmert.comIntroduction to Associated Legendre Polynomials (Legendre.jl Series, Part I)\nCalculating Legendre Polynomials (Legendre.jl Series, Part II)\nPre-normalizing Legendre Polynomials (Legendre.jl Series, Part III)\nMaintaining numerical accuracy in the Legendre recurrences (Legendre.jl Series, Part IV)\nIntroducing Legendre.jl (Legendre.jl Series, Part V)\nNumerical Accuracy of the Spherical Harmonic Recurrence Coefficient (Legendre.jl Series Addendum)\nNotes on Calculating the Spherical Harmonics\nMore Notes on Calculating the Spherical Harmonics: Analysis of maps to harmonic coefficients","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[Daley78]: Roger Daley & Yvon Bourassa (1978) Rhomboidal versus triangular spherical harmonic truncation: Some verification statistics, Atmosphere-Ocean, 16:2, 187-196, DOI: 10.1080/07055900.1978.9649026","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[Randall2021]: David Randall, 2021. An Introduction to Numerical Modeling of the Atmosphere, Chapter 22.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[Durran2010]: Dale Durran, 2010. Numerical Methods for Fluid Dynamics, Springer. In particular section 6.2, 6.4.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[GFDL]: Geophysical Fluid Dynamics Laboratory, The barotropic vorticity equation.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[FFT]: Depending on the implementation of the Fast Fourier Transform (Cooley-Tukey algorithm, or or the Bluestein algorithm) easily Fourier-transformable can mean different things: Vectors of the length n that is a power of two, i.e. n = 2^i is certainly easily Fourier-transformable, but for most FFT implementations so are n = 2^i3^j5^k with i j k some positive integers. In fact, FFTW uses O(n log n) algorithms even for prime sizes.","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[Bourke72]: Bourke, W. An Efficient, One-Level, Primitive-Equation Spectral Model. Mon. Wea. Rev. 100, 683–689 (1972). doi:10.1175/1520-0493(1972)100<0683:AEOPSM>2.3.CO;2","category":"page"},{"location":"spectral_transform/","page":"Spherical Harmonic Transform","title":"Spherical Harmonic Transform","text":"[Orszag70]: Orszag, S. A., 1970: Transform Method for the Calculation of Vector-Coupled Sums: Application to the Spectral Form of the Vorticity Equation. J. Atmos. Sci., 27, 890–895, 10.1175/1520-0469(1970)027<0890:TMFTCO>2.0.CO;2. ","category":"page"},{"location":"callbacks/#Callbacks","page":"Callbacks","title":"Callbacks","text":"","category":"section"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"SpeedyWeather.jl implements a callback system to let users include a flexible piece of code into the time stepping. You can think about the main time loop calling back to check whether anything else should be done before continuing with the next time step. The callback system here is called after the time step only (plus one call at initialize! and one at finalize!), we currently do not implement other callsites.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Callbacks are mainly introduced for diagnostic purposes, meaning that they do not influence the simulation, and access the prognostic variables and the model components in a read-only fashion. However, a callback is not strictly prevented from changing prognostic or diagnostic variables or the model. For example, you may define a callback that changes the orography during the simulation. In general, one has to keep the general order of executions during a time step in mind (valid for all models)","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"set tendencies to zero\ncompute parameterizations, forcing, or drag terms. Accumulate tendencies.\ncompute dynamics, accumulate tendencies.\ntime stepping\noutput\ncallbacks","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"This means that, at the current callsite, a callback can read the tendencies but writing into it would be overwritten by the zeroing of the tendencies in 1. anyway. At the moment, if a callback wants to implement an additional tendency then it currently should be implemented as a parameterization, forcing or drag term. ","category":"page"},{"location":"callbacks/#Defining-a-callback","page":"Callbacks","title":"Defining a callback","text":"","category":"section"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"You can (and are encouraged!) to write your own callbacks to diagnose SpeedyWeather simulations. Let us implement a StormChaser callback, recording the highest surface wind speed on every time step, that we want to use to illustrate how a callback needs to be defined.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Every custom callback needs to be defined as a (mutable) struct, subtype of AbstractCallback, i.e. struct or mutable struct CustomCallback <: SpeedyWeather.AbstractCallback. In our case, this is","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"using SpeedyWeather\n\nBase.@kwdef mutable struct StormChaser{NF} <: SpeedyWeather.AbstractCallback\n timestep_counter::Int = 0\n maximum_surface_wind_speed::Vector{NF} = [0]\nend\n\n# Generator function\nStormChaser(SG::SpectralGrid) = StormChaser{SG.NF}()","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"We decide to have a field timestep_counter in the callback that allows us to track the number of times the callback was called to create a time series of our highest surface wind speeds. The actual maximum_surface_wind_speed is then a vector of a given type NF (= number format), which is where we'll write into. Both are initialised with zeros. We also add a generator function, similar as to many other components in SpeedyWeather that just pulls the number format from the SpectralGrid object.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Now every callback needs to extend three methods","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"initialize!, called once before the main time loop starts\ncallback!, called after every time step\nfinalize!, called once after the last time step","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"And we'll go through them one by one.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"function SpeedyWeather.initialize!(\n callback::StormChaser,\n progn::PrognosticVariables,\n diagn::DiagnosticVariables,\n model::AbstractModel,\n)\n # allocate recorder: number of time steps (incl initial conditions) in simulation \n callback.maximum_surface_wind_speed = zeros(progn.clock.n_timesteps + 1)\n \n # where surface (=lowermost model layer) u, v on the grid are stored\n u_grid = diagn.grid.u_grid[:, diagn.nlayers]\n v_grid = diagn.grid.u_grid[:, diagn.nlayers]\n\n # maximum wind speed of initial conditions\n callback.maximum_surface_wind_speed[1] = max_2norm(u_grid, v_grid)\n \n # (re)set counter to 1\n callback.timestep_counter = 1\nend","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"The initialize! function has to be extended for the new callback ::StormChaser as first argument, then followed by prognostic and diagnostic variables and model. For correct multiple dispatch it is important to restrict the first argument to the new StormChaser type (to not call another callback instead), but the other type declarations are for clarity only. initialize!(::AbstractCallback, args...) is called once just before the main time loop, meaning after the initial conditions are set and after all other components are initialized. We replace the vector inside our storm chaser with a vector of the correct length so that we have a \"recorder\" allocated, a vector that can store the maximum surface wind speed on every time step. We then also compute that maximum for the initial conditions and set the time step counter to 1. We define the max_2norm function as follows","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"\"\"\"Maximum of the 2-norm of elements across two arrays.\"\"\"\nfunction max_2norm(u::AbstractArray{T}, v::AbstractArray{T}) where T\n max_norm = zero(T) # = u² + v²\n for ij in eachindex(u, v)\n # find largest wind speed squared\n max_norm = max(max_norm, u[ij]^2 + v[ij]^2)\n end\n return sqrt(max_norm) # take sqrt only once\nend","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Note that this function is defined in the scope Main and not inside SpeedyWeather, this is absolutely possible due to Julia's scope of variables which will use max_2norm from Main scope if it doesn't exist in the global scope inside the SpeedyWeather module scope. Then we need to extend the callback! function as follows","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"function SpeedyWeather.callback!(\n callback::StormChaser,\n progn::PrognosticVariables,\n diagn::DiagnosticVariables,\n model::AbstractModel,\n)\n\n # increase counter\n callback.timestep_counter += 1 \n i = callback.timestep_counter\n\n # where surface (=lowermost model layer) u, v on the grid are stored\n u_grid = diagn.grid.u_grid[:, diagn.nlayers]\n v_grid = diagn.grid.u_grid[:, diagn.nlayers]\n\n # maximum wind speed at current time step\n callback.maximum_surface_wind_speed[i] = max_2norm(u_grid, v_grid)\nend","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"The function signature for callback! is the same as for initialize!. You may access anything from progn, diagn or model, although for a purely diagnostic callback this should be read-only. While you could change other model components like the land sea mask in model.land_sea_mask or orography etc. then you interfere with the simulation which is more advanced and will be discussed in Intrusive callbacks below.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Lastly, we extend the finalize! function which is called once after the last time step. This could be used, for example, to save the maximum_surface_wind_speed vector to file or in case you want to find the highest wind speed across all time steps. But in many cases you may not need to do anything, in which case you just just let it return nothing.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"SpeedyWeather.finalize!(::StormChaser, args...) = nothing","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"note: Always extend `initialize!`, `callback!` and `finalize!`\nFor a custom callback you need to extend all three, initialize!, callback! and finalize!, even if your callback doesn't need it. Just return nothing in that case. Otherwise a MethodError will occur. While we could have defined all callbacks by default to do nothing on each of these, this may give you the false impression that your callback is already defined correctly, although it's not.","category":"page"},{"location":"callbacks/#Adding-a-callback","page":"Callbacks","title":"Adding a callback","text":"","category":"section"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Every model has a field callbacks::Dict{Symbol, AbstractCallback} such that the callbacks keyword can be used to create a model with a dictionary of callbacks. Callbacks are identified with a Symbol key inside such a dictionary. We have a convenient CallbackDict generator function which can be used like Dict but the key-value pairs have to be of type Symbol-AbstractCallback. Let us illustrate this with the dummy callback NoCallback (which is a callback that returns nothing on initialize!, callback! and finalize!)","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"callbacks = CallbackDict() # empty dictionary\ncallbacks = CallbackDict(:my_callback => NoCallback()) # key => callback","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"If you don't provide a key a random key will be assigned","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"callbacks = CallbackDict(NoCallback())","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"and you can add (or delete) additional callbacks","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"add!(callbacks, NoCallback()) # this will also pick a random key\nadd!(callbacks, :my_callback => NoCallback()) # use key :my_callback\ndelete!(callbacks, :my_callback) # remove by key\ncallbacks","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"And you can chain them too","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"add!(callbacks, NoCallback(), NoCallback()) # random keys\nadd!(callbacks, :key1 => NoCallback(), :key2 => NoCallback()) # keys provided","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Meaning that callbacks can be added before and after model construction","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"spectral_grid = SpectralGrid()\ncallbacks = CallbackDict(:callback_added_before => NoCallback())\nmodel = PrimitiveWetModel(; spectral_grid, callbacks)\nadd!(model.callbacks, :callback_added_afterwards => NoCallback())\nadd!(model, :callback_added_afterwards2 => NoCallback())","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Note how the first argument can be model.callbacks as outlined in the sections above because this is the callbacks dictionary, but also simply model, which will add the callback to model.callbacks. It's equivalent. Let us add two more meaningful callbacks","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"storm_chaser = StormChaser(spectral_grid)\nrecord_surface_temperature = GlobalSurfaceTemperatureCallback(spectral_grid)\nadd!(model.callbacks, :storm_chaser => storm_chaser)\nadd!(model.callbacks, :temperature => record_surface_temperature)","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"which means that now in the calls to callback! first the two dummy NoCallbacks are called and then our storm chaser callback and then the GlobalSurfaceTemperatureCallback which records the global mean surface temperature on every time step. From normal NetCDF output the information these callbacks analyse would not be available, only at the frequency of the model output, which for every time step would create way more data and considerably slow down the simulation. Let's run the simulation and check the callbacks","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"simulation = initialize!(model)\nrun!(simulation, period=Day(3))\nv = model.callbacks[:storm_chaser].maximum_surface_wind_speed\nmaximum(v) # highest surface wind speeds in simulation [m/s]","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Cool, our StormChaser callback with the key :storm_chaser has been recording maximum surface wind speeds in [m/s]. And the :temperature callback a time series of global mean surface temperatures in Kelvin on every time step while the model ran for 3 days.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"model.callbacks[:temperature].temp","category":"page"},{"location":"callbacks/#intrusive_callbacks","page":"Callbacks","title":"Intrusive callbacks","text":"","category":"section"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"In the sections above, callbacks were introduced as a tool to define custom diagnostics or simulation output. This is the simpler and recommended way of using them but nothing stops you from defining a callback that is intrusive, meaning that it can alter the prognostic or diagnostic variables or the model.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Changing any components of the model, e.g. boundary conditions like orography or the land-sea mask through a callback is possible although one should notice that this only comes into effect on the next time step given the execution order mentioned above. One could for example run a simulation for a certain period and then start moving continents around. Note that for physical consistency this should be reflected in the orography, land-sea mask, as well as in the available sea and land-surface temperatures, but one is free to do this only partially too. Another example would be to switch on/off certain model components over time. If these components are implemented as mutable struct then one could define a callback that weakens their respective strength parameter over time.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"As an example of a callback that changes the model components see","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Millenium flood: Time-dependent land-sea mask","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Changing the diagnostic variables, however, will not have any effect. All of them are treated as work arrays, meaning that their state is completely overwritten on every time step. Changing the prognostic variables in spectral space directly is not advised though possible because this can easily lead to stability issues. It is generally easier to implement something like this as a parameterization, forcing or drag term (which can also be made time-dependent).","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Overall, callbacks give the user a wide range of possibilities to diagnose the simulation while running or to interfere with a simulation. We therefore encourage users to use callbacks as widely as possible, but if you run into any issues please open an issue in the repository and explain what you'd like to achieve and which errors you are facing. We are happy to help.","category":"page"},{"location":"callbacks/#Schedules","page":"Callbacks","title":"Schedules","text":"","category":"section"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"For convenience, SpeedyWeather.jl implements a Schedule which helps to schedule when callbacks are called. Because in many situations you don't want to call them on every time step but only periodically, say once a day, or only on specific dates and times, e.g. Jan 1 at noon. Several examples how to create schedules","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"using SpeedyWeather\n\n# execute on timestep at or after Jan 2 2000\nevent_schedule = Schedule(DateTime(2000,1,2)) \n\n# several events scheduled\nevents = (DateTime(2000,1,3), DateTime(2000,1,5,12))\nseveral_events_schedule = Schedule(events...)\n\n# provided as Vector{DateTime} with times= keyword\nalways_at_noon = [DateTime(2000,1,i,12) for i in 1:10]\nnoon_schedule = Schedule(times=always_at_noon)\n\n# or using every= for periodic execution, here once a day\nperiodic_schedule = Schedule(every=Day(1))","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"A Schedule has 5 fields, see Schedule. every is an option to create a periodic schedule to execute every time that indicated period has passed. steps and counter will let you know how many callback execution steps there are and count them up. times is a Vector{DateTime} containing scheduled events. schedule is the actual schedule inside a Schedule, implemented as BitVector indicating whether to execute on a given time step (true) or not (false).","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Let's show how to use a Schedule inside a callback","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"struct MyScheduledCallback <: SpeedyWeather.AbstractCallback\n schedule::Schedule\n # add other fields here that you need\nend\n\nfunction SpeedyWeather.initialize!(\n callback::MyScheduledCallback,\n progn::PrognosticVariables,\n args...\n)\n # when initializing a scheduled callback also initialize its schedule!\n initialize!(callback.schedule, progn.clock)\n\n # initialize other things in your callback here\nend\n\nfunction SpeedyWeather.callback!(\n callback::MyScheduledCallback,\n progn::PrognosticVariables,\n diagn::DiagnosticVariables,\n model::AbstractModel,\n)\n # scheduled callbacks start with this line to execute only when scheduled!\n # else escape immediately\n isscheduled(callback.schedule, progn.clock) || return nothing\n\n # Just print the North Pole surface temperature to screen\n (;time) = progn.clock\n temp_at_north_pole = diagn.grid.temp_grid[1,end]\n\n @info \"North pole has a temperature of $temp_at_north_pole on $time.\"\nend\n\n# nothing needs to be done after simulation is finished\nSpeedyWeather.finalize!(::MyScheduledCallback, args...) = nothing","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"So in summary","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"add a field schedule::Schedule to your callback\nadd the line initialize!(callback.schedule, progn.clock) when initializing your callback\nstart your callback! method with isscheduled(callback.schedule, progn.clock) || return nothing to execute only when scheduled","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"A Schedule is a field inside a callback as this allows you the set the callbacks desired schedule when creating it. In the example above we can create our callback that is supposed to print the North Pole's temperature like so","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"north_pole_temp_at_noon_jan9 = MyScheduledCallback(Schedule(DateTime(2000,1,9,12)))","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"The default for every is typemax(Int) indicating \"never\". This just means that there is no periodically reoccuring schedule, only schedule.times would include some times for events that are scheduled. Now let's create a primitive equation model with that callback","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"spectral_grid = SpectralGrid(trunc=31, nlayers=5)\nmodel = PrimitiveWetModel(;spectral_grid)\nadd!(model.callbacks, north_pole_temp_at_noon_jan9)\n\n# start simulation 7 days earlier\nsimulation = initialize!(model, time = DateTime(2000,1,2,12))\nrun!(simulation, period=Day(10))\nnothing # hide","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"So the callback gives us the temperature at the North Pole exactly when scheduled. We could have also stored this temperature, or conditionally changed parameters inside the model. There are plenty of ways how to use the scheduling, either by event, or in contrast, we could also schedule for once a day. As illustrated in the following","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"north_pole_temp_daily = MyScheduledCallback(Schedule(every=Day(1)))\nadd!(model.callbacks, north_pole_temp_daily)\n\n# resume simulation, start time is now 2000-1-12 noon\nrun!(simulation, period=Day(5))\nnothing # hide","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Note that the previous callback is still part of the model, we haven't deleted it with delete!. But because it's scheduled for a specific time that's in the past now that we resume the simulation it's schedule is empty (which is thrown as a warning). However, our new callback, scheduled daily, is active and prints daily at noon, because the simulation start time was noon.","category":"page"},{"location":"callbacks/#Scheduling-logic","page":"Callbacks","title":"Scheduling logic","text":"","category":"section"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"An event Schedule (created with DateTime object(s)) for callbacks, executes on or after the specified times. For two consecutive time steps i, i+1, an event is scheduled at i+1 when it occurs in (ii+1. So a simulation with timestep i on Jan-1 at 1am, and i+1 at 2am, will execute a callback scheduled for 1am at i but scheduled for 1am and 1s (=01:00:01 on a 24H clock) at 2am. Because callbacks are always executed after a timestep this also means that a simulation starting at midnight with a callback scheduled for midnight will not execute this callback as it is outside of the (i i+1 range. You'd need to include this execution into the initialization. If several events inside the Schedule fall into the same time step (in the example above, 1am and 1s and 1am 30min) the execution will not happen twice. Think of a scheduled callback as a binary \"should the callback be executed now or not?\". Which is in fact how it's implemented, as a BitVector of the length of the number of time steps. If the bit at a given timestep is true, execute, otherwise not.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"A periodic Schedule (created with every = Hour(2) or similar) will execute on the timestep after that period (here 2 hours) has passed. If a simulation starts at midnight with one hour time steps then execution would take place after the timestep from 1am to 2am because that's when the clock switches to 2am which is 2 hours after the start of the simulation. Note that therefore the initial timestep is not included, however, the last time step would be if the period is a multiple of the scheduling period. If the first timestep should be included (e.g. you want to do something with the initial conditions) then you'll need to include that into the initialization of the callback.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Periodic schedules which do not match the simulation time step will be adjusted by rounding. Example, if you want a schedule which executes every hour but your simulation time step is 25min then it will be adjusted to execute every 2nd time step, meaning every 50min and not 1 hour. However, an info will be thrown if that is the case","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"odd_schedule = MyScheduledCallback(Schedule(every = Minute(70)))\nadd!(model.callbacks, odd_schedule)\n\n# resume simulation for 4 hours\nrun!(simulation, period=Hour(4))\nnothing # hide","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Now we get two empty schedules, one from callback that's supposed to execute on Jan 9 noon (this time has passed in our simulation) and one from the daily callback (we're not simulating for a day). You could just delete! those callbacks. You can see that while we wanted our odd_schedule to execute every 70min, it has to adjust it to every 60min to match the simulation time step of 30min.","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"After the model initialization you can always check the simulation time step from model.time_stepping ","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"model.time_stepping","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Or converted into minutes (the time step internally is at millisecond accuracy)","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"Minute(model.time_stepping.Δt_millisec)","category":"page"},{"location":"callbacks/","page":"Callbacks","title":"Callbacks","text":"which illustrates why the adjustment of our callback frequency was necessary.","category":"page"},{"location":"ringgrids/#RingGrids","page":"RingGrids","title":"RingGrids","text":"","category":"section"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"RingGrids is a submodule that has been developed for SpeedyWeather.jl which is technically independent (SpeedyWeather.jl however imports it and so does SpeedyTransforms) and can also be used without running simulations. It is just not put into its own respective repository.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"RingGrids defines several iso-latitude grids, which are mathematically described in the section on Grids. In brief, they include the regular latitude-longitude grids (here called FullClenshawGrid) as well as grids which latitudes are shifted to the Gaussian latitudes and reduced grids, meaning that they have a decreasing number of longitudinal points towards the poles to be more equal-area than full grids.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"RingGrids defines and exports the following grids:","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"full grids: FullClenshawGrid, FullGaussianGrid, FullHEALPix, and FullOctaHEALPix\nreduced grids: OctahedralGaussianGrid, OctahedralClenshawGrid, OctaHEALPixGrid and HEALPixGrid","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"The following explanation of how to use these can be mostly applied to any of them, however, there are certain functions that are not defined, e.g. the full grids can be trivially converted to a Matrix (i.e. they are rectangular grids) but not the OctahedralGaussianGrid.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"note: What is a ring?\nWe use the term ring, short for iso-latitude ring, to refer to a sequence of grid points that all share the same latitude. A latitude-longitude grid is a ring grid, as it organises its grid-points into rings. However, other grids, like the cubed-sphere are not based on iso-latitude rings. SpeedyWeather.jl only works with ring grids because its a requirement for the Spherical Harmonic Transform to be efficient. See Grids.","category":"page"},{"location":"ringgrids/#Creating-data-on-a-RingGrid","page":"RingGrids","title":"Creating data on a RingGrid","text":"","category":"section"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"Every grid in RingGrids has a grid.data field, which is a vector containing the data on the grid. The grid points are unravelled west to east then north to south, meaning that it starts at 90˚N and 0˚E then walks eastward for 360˚ before jumping on the next latitude ring further south, this way circling around the sphere till reaching the south pole. This may also be called ring order.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"Data in a Matrix which follows this ring order can be put on a FullGaussianGrid like so","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"using SpeedyWeather.RingGrids\nmap = randn(Float32, 8, 4)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"grid = FullGaussianGrid(map, input_as=Matrix)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"Note that input_as=Matrix is necessary as, RingGrids have a flattened horizontal dimension into a vector. To distinguish the 2nd horizontal dimension from a possible vertical dimension the keyword argument here is required.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"A full Gaussian grid has always 2N x N grid points, but a FullClenshawGrid has 2N x N-1, if those dimensions don't match, the creation will throw an error. To reobtain the data from a grid, you can access its data field which returns a normal Vector","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"grid.data","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"Which can be reshaped to reobtain map from above. Alternatively you can Matrix(grid) to do this in one step","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"map == Matrix(FullGaussianGrid(map))","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"You can also use zeros, ones, rand, randn to create a grid, whereby nlat_half, i.e. the number of latitude rings on one hemisphere, Equator included, is used as a resolution parameter and here as a second argument.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"nlat_half = 4\ngrid = randn(OctahedralGaussianGrid{Float16}, nlat_half)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"and any element type T can be used for OctahedralGaussianGrid{T} and similar for other grid types.","category":"page"},{"location":"ringgrids/#Visualising-RingGrid-data","page":"RingGrids","title":"Visualising RingGrid data","text":"","category":"section"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"As only the full grids can be reshaped into a matrix, the underlying data structure of any AbstractGrid is a vector. As shown in the examples above, one can therefore inspect the data as if it was a vector. But as that data has, through its <:AbstractGrid type, all the geometric information available to plot it on a map, RingGrids also exports plot function, based on UnicodePlots' heatmap.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"nlat_half = 24\ngrid = randn(OctahedralGaussianGrid, nlat_half)\nRingGrids.plot(grid)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"(Note that to skip the RingGrids. in the last line you can do import SpeedyWeather.RingGrids: plot, import SpeedyWeather: plot or simply using SpeedyWeather. It's just that LowerTriangularMatrices also defines plot which otherwise causes naming conflicts.)","category":"page"},{"location":"ringgrids/#Indexing-RingGrids","page":"RingGrids","title":"Indexing RingGrids","text":"","category":"section"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"All RingGrids have a single index ij which follows the ring order. While this is obviously not super exciting here are some examples how to make better use of the information that the data sits on a grid.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"We obtain the latitudes of the rings of a grid by calling get_latd (get_lond is only defined for full grids, or use get_latdlonds for latitudes, longitudes per grid point not per ring)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"grid = randn(OctahedralClenshawGrid, 5)\nlatd = get_latd(grid)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"Now we could calculate Coriolis and add it on the grid as follows","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"rotation = 7.29e-5 # angular frequency of Earth's rotation [rad/s]\ncoriolis = 2rotation*sind.(latd) # vector of coriolis parameters per latitude ring\n\nrings = eachring(grid)\nfor (j, ring) in enumerate(rings)\n f = coriolis[j]\n for ij in ring\n grid[ij] += f\n end\nend","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"eachring creates a vector of UnitRange indices, such that we can loop over the ring index j (j=1 being closest to the North pole) pull the coriolis parameter at that latitude and then loop over all in-ring indices i (changing longitudes) to do something on the grid. Something similar can be done to scale/unscale with the cosine of latitude for example. We can always loop over all grid-points like so","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"for ij in eachgridpoint(grid)\n grid[ij]\nend","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"or use eachindex instead.","category":"page"},{"location":"ringgrids/#Interpolation-on-RingGrids","page":"RingGrids","title":"Interpolation on RingGrids","text":"","category":"section"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"In most cases we will want to use RingGrids so that our data directly comes with the geometric information of where the grid-point is one the sphere. We have seen how to use get_latd, get_lond, ... for that above. This information generally can also be used to interpolate our data from grid to another or to request an interpolated value on some coordinates. Using our data on grid which is an OctahedralGaussianGrid from above we can use the interpolate function to get it onto a FullGaussianGrid (or any other grid for purpose)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"grid = randn(OctahedralGaussianGrid{Float32}, 4)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"interpolate(FullGaussianGrid, grid)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"By default this will linearly interpolate (it's an Anvil interpolator, see below) onto a grid with the same nlat_half, but we can also coarse-grain or fine-grain by specifying nlat_half directly as 2nd argument","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"interpolate(FullGaussianGrid, 6, grid)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"So we got from an 8-ring OctahedralGaussianGrid{Float16} to a 12-ring FullGaussianGrid{Float64}, so it did a conversion from Float16 to Float64 on the fly too, because the default precision is Float64 unless specified. interpolate(FullGaussianGrid{Float16}, 6, grid) would have interpolated onto a grid with element type Float16.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"One can also interpolate onto a given coordinate ˚N, ˚E like so","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"interpolate(30.0, 10.0, grid)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"we interpolated the data from grid onto 30˚N, 10˚E. To do this simultaneously for many coordinates they can be packed into a vector too","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"interpolate([30.0, 40.0, 50.0], [10.0, 10.0, 10.0], grid)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"which returns the data on grid at 30˚N, 40˚N, 50˚N, and 10˚E respectively. Note how the interpolation here retains the element type of grid.","category":"page"},{"location":"ringgrids/#Performance-for-RingGrid-interpolation","page":"RingGrids","title":"Performance for RingGrid interpolation","text":"","category":"section"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"Every time an interpolation like interpolate(30.0, 10.0, grid) is called, several things happen, which are important to understand to know how to get the fastest interpolation out of this module in a given situation. Under the hood an interpolation takes three arguments","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"output vector\ninput grid\ninterpolator","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"The output vector is just an array into which the interpolated data is written, providing this prevents unnecessary allocation of memory in case the destination array of the interpolation already exists. The input grid contains the data which is subject to interpolation, it must come on a ring grid, however, its coordinate information is actually already in the interpolator. The interpolator knows about the geometry of the grid the data is coming on and the coordinates it is supposed to interpolate onto. It has therefore precalculated the indices that are needed to access the right data on the input grid and the weights it needs to apply in the actual interpolation operation. The only thing it does not know is the actual data values of that grid. So in the case you want to interpolate from grid A to grid B many times, you can just reuse the same interpolator. If you want to change the coordinates of the output grid but its total number of points remain constants then you can update the locator inside the interpolator and only else you will need to create a new interpolator. Let's look at this in practice. Say we have two grids an want to interpolate between them","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"grid_in = rand(HEALPixGrid, 4)\ngrid_out = zeros(FullClenshawGrid, 6)\ninterp = RingGrids.interpolator(grid_out, grid_in)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"Now we have created an interpolator interp which knows about the geometry where to interpolate from and the coordinates there to interpolate to. It is also initialized, meaning it has precomputed the indices to of grid_in that are supposed to be used. It just does not know about the data of grid_in (and neither of grid_out which will be overwritten anyway). We can now do","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"interpolate!(grid_out, grid_in, interp)\ngrid_out","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"which is identical to interpolate(grid_out, grid_in) but you can reuse interp for other data. The interpolation can also handle various element types (the interpolator interp does not have to be updated for this either)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"grid_out = zeros(FullClenshawGrid{Float16}, 6);\ninterpolate!(grid_out, grid_in, interp)\ngrid_out","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"and we have converted data from a HEALPixGrid{Float64} (Float64 is always default if not specified) to a FullClenshawGrid{Float16} including the type conversion Float64-Float16 on the fly. Technically there are three data types and their combinations possible: The input data will come with a type, the output array has an element type and the interpolator has precomputed weights with a given type. Say we want to go from Float16 data on an OctahedralGaussianGrid to Float16 on a FullClenshawGrid but using Float32 precision for the interpolation itself, we would do this by","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"grid_in = randn(OctahedralGaussianGrid{Float16}, 24)\ngrid_out = zeros(FullClenshawGrid{Float16}, 24)\ninterp = RingGrids.interpolator(Float32, grid_out, grid_in)\ninterpolate!(grid_out, grid_in, interp)\ngrid_out","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"As a last example we want to illustrate a situation where we would always want to interpolate onto 10 coordinates, but their locations may change. In order to avoid recreating an interpolator object we would do (AnvilInterpolator is described in Anvil interpolator)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"npoints = 10 # number of coordinates to interpolate onto\ninterp = AnvilInterpolator(Float32, HEALPixGrid, 24, npoints)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"with the first argument being the number format used during interpolation, then the input grid type, its resolution in terms of nlat_half and then the number of points to interpolate onto. However, interp is not yet initialized as it does not know about the destination coordinates yet. Let's define them, but note that we already decided there's only 10 of them above.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"latds = collect(0.0:5.0:45.0)\nlonds = collect(-10.0:2.0:8.0)\nnothing # hide","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"now we can update the locator inside our interpolator as follows","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"RingGrids.update_locator!(interp, latds, londs)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"With data matching the input from above, a nlat_half=24 HEALPixGrid, and allocate 10-element output vector","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"output_vec = zeros(10)\ngrid_input = rand(HEALPixGrid, 24)\nnothing # hide","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"we can use the interpolator as follows","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"interpolate!(output_vec, grid_input, interp)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"which is the approximately the same as doing it directly without creating an interpolator first and updating its locator","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"interpolate(latds, londs, grid_input)","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"but allows for a reuse of the interpolator. Note that the two output arrays are not exactly identical because we manually set our interpolator interp to use Float32 for the interpolation whereas the default is Float64.","category":"page"},{"location":"ringgrids/#Anvil-interpolator","page":"RingGrids","title":"Anvil interpolator","text":"","category":"section"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"Currently the only interpolator implemented is a 4-point bilinear interpolator, which schematically works as follows. Anvil interpolation is the bilinear average of a, b, c, d which are values at grid points in an anvil-shaped configuration at location x, which is denoted by Δab, Δcd, Δy, the fraction of distances between a-b, c-d, and ab-cd, respectively. Note that a, c and b, d do not necessarily share the same longitude/x-coordinate.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":" 0..............1 # fraction of distance Δab between a, b\n |< Δab >|\n\n0^ a -------- o - b # anvil-shaped average of a, b, c, d at location x\n.Δy |\n. |\n.v x \n. |\n1 c ------ o ---- d\n\n |< Δcd >|\n 0...............1 # fraction of distance Δcd between c, d\n\n^ fraction of distance Δy between a-b and c-d.","category":"page"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"This interpolation is chosen as by definition of the ring grids, a and b share the same latitude, so do c and d, but the longitudes can be different for all four, a, b, c, d.","category":"page"},{"location":"ringgrids/#Function-index","page":"RingGrids","title":"Function index","text":"","category":"section"},{"location":"ringgrids/","page":"RingGrids","title":"RingGrids","text":"Modules = [SpeedyWeather.RingGrids]","category":"page"},{"location":"ringgrids/#Core.Type-Tuple{AbstractArray}","page":"RingGrids","title":"Core.Type","text":"() Initialize an instance of the grid from an Array. For keyword argument input_as=Vector (default) the leading dimension is interpreted as a flat vector of all horizontal entries in one layer. For input_as==Matrx the first two leading dimensions are interpreted as longitute and latitude. This is only possible for full grids that are a subtype of AbstractFullGridArray.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AbstractFullGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.AbstractFullGrid","text":"An AbstractFullGrid is a horizontal grid with a constant number of longitude points across latitude rings. Different latitudes can be used, Gaussian latitudes, equi-angle latitudes (also called Clenshaw from Clenshaw-Curtis quadrature), or others.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AbstractFullGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.AbstractFullGridArray","text":"Subtype of AbstractGridArray for all N-dimensional arrays of ring grids that have the same number of longitude points on every ring. As such these (horizontal) grids are representable as a matrix, with denser grid points towards the poles.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AbstractGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.AbstractGrid","text":"Abstract supertype for all ring grids, representing 2-dimensional data on the sphere unravelled into a Julia Vector. Subtype of AbstractGridArray with N=1 and ArrayType=Vector{T} of eltype T.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.AbstractGridArray","text":"Abstract supertype for all arrays of ring grids, representing N-dimensional data on the sphere in two dimensions (but unravelled into a vector in the first dimension, the actual \"ring grid\") plus additional N-1 dimensions for the vertical and/or time etc. Parameter T is the eltype of the underlying data, held as in the array type ArrayType (Julia's Array for CPU or others for GPU).\n\nRing grids have several consecuitive grid points share the same latitude (= a ring), grid points on a given ring are equidistant. Grid points are ordered 0 to 360˚E, starting around the north pole, ring by ring to the south pole. \n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AbstractInterpolator","page":"RingGrids","title":"SpeedyWeather.RingGrids.AbstractInterpolator","text":"abstract type AbstractInterpolator{NF, G} end\n\nSupertype for Interpolators. Every Interpolator <: AbstractInterpolator is expected to have two fields,\n\ngeometry, which describes the grid G to interpolate from\nlocator, which locates the indices on G and their weights to interpolate onto a new grid.\n\nNF is the number format used to calculate the interpolation, which can be different from the input data and/or the interpolated data on the new grid.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AbstractLocator","page":"RingGrids","title":"SpeedyWeather.RingGrids.AbstractLocator","text":"AbstractLocator{NF}\n\nSupertype of every Locator, which locates the indices on a grid to be used to perform an interpolation. E.g. AnvilLocator uses a 4-point stencil for every new coordinate to interpolate onto. Higher order stencils can be implemented by defining OtherLocator <: AbstractLocactor.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AbstractReducedGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.AbstractReducedGrid","text":"Horizontal abstract type for all AbstractReducedGridArray with N=1 (i.e. horizontal only) and ArrayType of Vector{T} with element type T.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AbstractReducedGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.AbstractReducedGridArray","text":"Subtype of AbstractGridArray for arrays of rings grids that have a reduced number of longitude points towards the poles, i.e. they are not \"full\", see AbstractFullGridArray. Data on these grids cannot be represented as matrix and has to be unravelled into a vector, ordered 0 to 360˚E then north to south, ring by ring. Examples for reduced grids are the octahedral Gaussian or Clenshaw grids, or the HEALPix grid.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AbstractSphericalDistance","page":"RingGrids","title":"SpeedyWeather.RingGrids.AbstractSphericalDistance","text":"abstract type AbstractSphericalDistance end\n\nSuper type of formulas to calculate the spherical distance or great-circle distance. To define a NewFormula, define struct NewFormula <: AbstractSphericalDistance end and the actual calculation as a functor\n\nfunction NewFormula(lonlat1::Tuple, lonlat2::Tuple; radius=DEFAULT_RADIUS, kwargs...)\n\nassuming inputs in degrees and returning the distance in meters (or radians for radius=1). Then use the general interface spherical_distance(NewFormula, args...; kwargs...)\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AnvilLocator","page":"RingGrids","title":"SpeedyWeather.RingGrids.AnvilLocator","text":"AnvilLocator{NF<:AbstractFloat} <: AbtractLocator\n\nContains arrays that locates grid points of a given field to be uses in an interpolation and their weights. This Locator is a 4-point average in an anvil-shaped grid-point arrangement between two latitude rings.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.AnvilLocator-Union{Tuple{Integer}, Tuple{NF}} where NF<:AbstractFloat","page":"RingGrids","title":"SpeedyWeather.RingGrids.AnvilLocator","text":"L = AnvilLocator( ::Type{NF}, # number format used for the interpolation\n npoints::Integer # number of points to interpolate onto\n ) where {NF<:AbstractFloat}\n\nZero generator function for the 4-point average AnvilLocator. Use update_locator! to update the grid indices used for interpolation and their weights. The number format NF is the format used for the calculations within the interpolation, the input data and/or output data formats may differ.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.FullClenshawArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.FullClenshawArray","text":"A FullClenshawArray is an array of full grid, subtyping AbstractFullGridArray, that use equidistant latitudes for each ring (a regular lon-lat grid). These require the Clenshaw-Curtis quadrature in the spectral transform, hence the name. One ring is on the equator, total number of rings is odd, no rings on the north or south pole. First dimension of the underlying N-dimensional data represents the horizontal dimension, in ring order (0 to 360˚E, then north to south), other dimensions are used for the vertical and/or time or other dimensions. The resolution parameter of the horizontal grid is nlat_half (number of latitude rings on one hemisphere, Equator included) and the ring indices are precomputed in rings. Fields are\n\ndata::AbstractArray\nnlat_half::Int64\nrings::Vector{UnitRange{Int64}}\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.FullClenshawGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.FullClenshawGrid","text":"A FullClenshawArray but constrained to N=1 dimensions (horizontal only) and data is a Vector{T}.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.FullGaussianArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.FullGaussianArray","text":"A FullGaussianArray is an array of full grids, subtyping AbstractFullGridArray, that use Gaussian latitudes for each ring. First dimension of the underlying N-dimensional data represents the horizontal dimension, in ring order (0 to 360˚E, then north to south), other dimensions are used for the vertical and/or time or other dimensions. The resolution parameter of the horizontal grid is nlat_half (number of latitude rings on one hemisphere, Equator included) and the ring indices are precomputed in rings. Fields are\n\ndata::AbstractArray: Data array, west to east, ring by ring, north to south.\nnlat_half::Int64: Number of latitudes on one hemisphere\nrings::Vector{UnitRange{Int64}}: Precomputed ring indices, ranging from first to last grid point on every ring.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.FullGaussianGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.FullGaussianGrid","text":"A FullGaussianArray but constrained to N=1 dimensions (horizontal only) and data is a Vector{T}.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.FullHEALPixArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.FullHEALPixArray","text":"A FullHEALPixArray is an array of full grids, subtyping AbstractFullGridArray, that use HEALPix latitudes for each ring. This type primarily equists to interpolate data from the (reduced) HEALPixGrid onto a full grid for output.\n\nFirst dimension of the underlying N-dimensional data represents the horizontal dimension, in ring order (0 to 360˚E, then north to south), other dimensions are used for the vertical and/or time or other dimensions. The resolution parameter of the horizontal grid is nlat_half (number of latitude rings on one hemisphere, Equator included) and the ring indices are precomputed in rings. Fields are\n\ndata::AbstractArray\nnlat_half::Int64\nrings::Vector{UnitRange{Int64}}\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.FullHEALPixGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.FullHEALPixGrid","text":"A FullHEALPixArray but constrained to N=1 dimensions (horizontal only) and data is a Vector{T}.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.FullOctaHEALPixArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.FullOctaHEALPixArray","text":"A FullOctaHEALPixArray is an array of full grids, subtyping AbstractFullGridArray that use OctaHEALPix latitudes for each ring. This type primarily equists to interpolate data from the (reduced) OctaHEALPixGrid onto a full grid for output.\n\nFirst dimension of the underlying N-dimensional data represents the horizontal dimension, in ring order (0 to 360˚E, then north to south), other dimensions are used for the vertical and/or time or other dimensions. The resolution parameter of the horizontal grid is nlat_half (number of latitude rings on one hemisphere, Equator included) and the ring indices are precomputed in rings. Fields are\n\ndata::AbstractArray\nnlat_half::Int64\nrings::Vector{UnitRange{Int64}}\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.FullOctaHEALPixGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.FullOctaHEALPixGrid","text":"A FullOctaHEALPixArray but constrained to N=1 dimensions (horizontal only) and data is a Vector{T}.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.GridGeometry","page":"RingGrids","title":"SpeedyWeather.RingGrids.GridGeometry","text":"GridGeometry{G<:AbstractGrid}\n\ncontains general precomputed arrays describing the grid of G.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.GridGeometry-Tuple{Type{<:AbstractGrid}, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.GridGeometry","text":"G = GridGeometry( Grid::Type{<:AbstractGrid},\n nlat_half::Integer)\n\nPrecomputed arrays describing the geometry of the Grid with resolution nlat_half. Contains latitudes and longitudes of grid points, their ring index j and their unravelled indices ij.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.HEALPixArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.HEALPixArray","text":"A HEALPixArray is an array of HEALPix grids, subtyping AbstractReducedGridArray. First dimension of the underlying N-dimensional data represents the horizontal dimension, in ring order (0 to 360˚E, then north to south), other dimensions are used for the vertical and/or time or other dimensions. The resolution parameter of the horizontal grid is nlat_half (number of latitude rings on one hemisphere, Equator included) which has to be even (non-fatal error thrown otherwise) which is less strict than the original HEALPix formulation (only powers of two for nside = nlat_half/2). Ring indices are precomputed in rings.\n\nA HEALPix grid has 12 faces, each nsidexnside grid points, each covering the same area of the sphere. They start with 4 longitude points on the northern-most ring, increase by 4 points per ring in the \"polar cap\" (the top half of the 4 northern-most faces) but have a constant number of longitude points in the equatorial belt. The southern hemisphere is symmetric to the northern, mirrored around the Equator. HEALPix grids have a ring on the Equator. For more details see Górski et al. 2005, DOI:10.1086/427976. \n\nrings are the precomputed ring indices, for nlat_half = 4 it is rings = [1:4, 5:12, 13:20, 21:28, 29:36, 37:44, 45:48]. So the first ring has indices 1:4 in the unravelled first dimension, etc. For efficient looping see eachring and eachgrid. Fields are\n\ndata::AbstractArray\nnlat_half::Int64\nrings::Vector{UnitRange{Int64}}\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.HEALPixGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.HEALPixGrid","text":"An HEALPixArray but constrained to N=1 dimensions (horizontal only) and data is a Vector{T}.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.Haversine-Tuple{Tuple, Tuple}","page":"RingGrids","title":"SpeedyWeather.RingGrids.Haversine","text":"Haversine(lonlat1::Tuple, lonlat2::Tuple; radius) -> Any\n\n\nHaversine formula calculating the great-circle or spherical distance (in meters) on the sphere between two tuples of longitude-latitude points in degrees ˚E, ˚N. Use keyword argument radius to change the radius of the sphere (default 6371e3 meters, Earth's radius), use radius=1 to return the central angle in radians or radius=360/2π to return degrees.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.OctaHEALPixArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.OctaHEALPixArray","text":"An OctaHEALPixArray is an array of OctaHEALPix grids, subtyping AbstractReducedGridArray. First dimension of the underlying N-dimensional data represents the horizontal dimension, in ring order (0 to 360˚E, then north to south), other dimensions are used for the vertical and/or time or other dimensions. The resolution parameter of the horizontal grid is nlat_half (number of latitude rings on one hemisphere, Equator included) and the ring indices are precomputed in rings.\n\nAn OctaHEALPix grid has 4 faces, each nlat_half x nlat_half in size, covering 90˚ in longitude, pole to pole. As part of the HEALPix family of grids, the grid points are equal area. They start with 4 longitude points on the northern-most ring, increase by 4 points per ring towards the Equator with one ring on the Equator before reducing the number of points again towards the south pole by 4 per ring. There is no equatorial belt for OctaHEALPix grids. The southern hemisphere is symmetric to the northern, mirrored around the Equator. OctaHEALPix grids have a ring on the Equator. For more details see Górski et al. 2005, DOI:10.1086/427976, the OctaHEALPix grid belongs to the family of HEALPix grids with Nθ = 1, Nφ = 4 but is not explicitly mentioned therein.\n\nrings are the precomputed ring indices, for nlat_half = 3 (in contrast to HEALPix this can be odd) it is rings = [1:4, 5:12, 13:24, 25:32, 33:36]. For efficient looping see eachring and eachgrid. Fields are\n\ndata::AbstractArray\nnlat_half::Int64\nrings::Vector{UnitRange{Int64}}\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.OctaHEALPixGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.OctaHEALPixGrid","text":"An OctaHEALPixArray but constrained to N=1 dimensions (horizontal only) and data is a Vector{T}.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.OctahedralClenshawArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.OctahedralClenshawArray","text":"An OctahedralClenshawArray is an array of octahedral grids, subtyping AbstractReducedGridArray, that use equidistant latitudes for each ring, the same as for FullClenshawArray. First dimension of the underlying N-dimensional data represents the horizontal dimension, in ring order (0 to 360˚E, then north to south), other dimensions are used for the vertical and/or time or other dimensions. The resolution parameter of the horizontal grid is nlat_half (number of latitude rings on one hemisphere, Equator included) and the ring indices are precomputed in rings.\n\nThese grids are called octahedral (same as for the OctahedralGaussianArray which only uses different latitudes) because after starting with 20 points on the first ring around the north pole (default) they increase the number of longitude points for each ring by 4, such that they can be conceptually thought of as lying on the 4 faces of an octahedron on each hemisphere. Hence, these grids have 20, 24, 28, ... longitude points for ring 1, 2, 3, ... Clenshaw grids have a ring on the Equator which has 16 + 4nlat_half longitude points before reducing the number of longitude points per ring by 4 towards the southern-most ring j = nlat. rings are the precomputed ring indices, the the example above rings = [1:20, 21:44, 45:72, ...]. For efficient looping see eachring and eachgrid. Fields are\n\ndata::AbstractArray\nnlat_half::Int64\nrings::Vector{UnitRange{Int64}}\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.OctahedralClenshawGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.OctahedralClenshawGrid","text":"An OctahedralClenshawArray but constrained to N=1 dimensions (horizontal only) and data is a Vector{T}.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.OctahedralGaussianArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.OctahedralGaussianArray","text":"An OctahedralGaussianArray is an array of octahedral grids, subtyping AbstractReducedGridArray, that use Gaussian latitudes for each ring. First dimension of the underlying N-dimensional data represents the horizontal dimension, in ring order (0 to 360˚E, then north to south), other dimensions are used for the vertical and/or time or other dimensions. The resolution parameter of the horizontal grid is nlat_half (number of latitude rings on one hemisphere, Equator included) and the ring indices are precomputed in rings.\n\nThese grids are called octahedral because after starting with 20 points on the first ring around the north pole (default) they increase the number of longitude points for each ring by 4, such that they can be conceptually thought of as lying on the 4 faces of an octahedron on each hemisphere. Hence, these grids have 20, 24, 28, ... longitude points for ring 1, 2, 3, ... There is no ring on the Equator and the two rings around it have 16 + 4nlat_half longitude points before reducing the number of longitude points per ring by 4 towards the southern-most ring j = nlat. rings are the precomputed ring indices, the the example above rings = [1:20, 21:44, 45:72, ...]. For efficient looping see eachring and eachgrid. Fields are\n\ndata::AbstractArray\nnlat_half::Int64\nrings::Vector{UnitRange{Int64}}\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#SpeedyWeather.RingGrids.OctahedralGaussianGrid","page":"RingGrids","title":"SpeedyWeather.RingGrids.OctahedralGaussianGrid","text":"An OctahedralGaussianArray but constrained to N=1 dimensions (horizontal only) and data is a Vector{T}.\n\n\n\n\n\n","category":"type"},{"location":"ringgrids/#Base.sizeof-Tuple{AbstractGridArray}","page":"RingGrids","title":"Base.sizeof","text":"sizeof(G::AbstractGridArray) -> Any\n\n\nSize of underlying data array plus precomputed ring indices.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.Matrix!-Tuple{AbstractMatrix, OctaHEALPixGrid}","page":"RingGrids","title":"SpeedyWeather.RingGrids.Matrix!","text":"Matrix!(M::AbstractMatrix,\n G::OctaHEALPixGrid;\n quadrant_rotation=(0, 1, 2, 3),\n matrix_quadrant=((2, 2), (1, 2), (1, 1), (2, 1)),\n )\n\nSorts the gridpoints in G into the matrix M without interpolation. Every quadrant of the grid G is rotated as specified in quadrant_rotation, 0 is no rotation, 1 is 90˚ clockwise, 2 is 180˚ etc. Grid quadrants are counted eastward starting from 0˚E. The grid quadrants are moved into the matrix quadrant (i, j) as specified. Defaults are equivalent to centered at 0˚E and a rotation such that the North Pole is at M's midpoint.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.Matrix!-Tuple{AbstractMatrix, OctahedralClenshawGrid}","page":"RingGrids","title":"SpeedyWeather.RingGrids.Matrix!","text":"Matrix!(\n M::AbstractMatrix,\n G::OctahedralClenshawGrid;\n kwargs...\n) -> AbstractMatrix\n\n\nSorts the gridpoints in G into the matrix M without interpolation. Every quadrant of the grid G is rotated as specified in quadrant_rotation, 0 is no rotation, 1 is 90˚ clockwise, 2 is 180˚ etc. Grid quadrants are counted eastward starting from 0˚E. The grid quadrants are moved into the matrix quadrant (i, j) as specified. Defaults are equivalent to centered at 0˚E and a rotation such that the North Pole is at M's midpoint.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.Matrix!-Union{Tuple{Vararg{Tuple{AbstractMatrix{T}, OctaHEALPixGrid}}}, Tuple{T}} where T","page":"RingGrids","title":"SpeedyWeather.RingGrids.Matrix!","text":"Matrix!(MGs::Tuple{AbstractMatrix{T}, OctaHEALPixGrid}...; kwargs...)\n\nLike Matrix!(::AbstractMatrix, ::OctaHEALPixGrid) but for simultaneous processing of tuples ((M1, G1), (M2, G2), ...) with matrices Mi and grids Gi. All matrices and grids have to be of the same size respectively.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.Matrix!-Union{Tuple{Vararg{Tuple{AbstractMatrix{T}, OctahedralClenshawGrid}}}, Tuple{T}} where T","page":"RingGrids","title":"SpeedyWeather.RingGrids.Matrix!","text":"Matrix!(\n MGs::Tuple{AbstractArray{T, 2}, OctahedralClenshawGrid}...;\n quadrant_rotation,\n matrix_quadrant\n) -> Union{Tuple, AbstractMatrix}\n\n\nLike Matrix!(::AbstractMatrix, ::OctahedralClenshawGrid) but for simultaneous processing of tuples ((M1, G1), (M2, G2), ...) with matrices Mi and grids Gi. All matrices and grids have to be of the same size respectively.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids._scale_lat!-Union{Tuple{T}, Tuple{AbstractGridArray{T, N, ArrayType} where {N, ArrayType<:AbstractArray{T, N}}, AbstractVector}} where T","page":"RingGrids","title":"SpeedyWeather.RingGrids._scale_lat!","text":"_scale_lat!(\n grid::AbstractGridArray{T, N, ArrayType} where {N, ArrayType<:AbstractArray{T, N}},\n v::AbstractVector\n) -> AbstractGridArray\n\n\nGeneric latitude scaling applied to A in-place with latitude-like vector v.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.anvil_average-NTuple{7, Any}","page":"RingGrids","title":"SpeedyWeather.RingGrids.anvil_average","text":"anvil_average(a, b, c, d, Δab, Δcd, Δy) -> Any\n\n\nThe bilinear average of a, b, c, d which are values at grid points in an anvil-shaped configuration at location x, which is denoted by Δab, Δcd, Δy, the fraction of distances between a-b, c-d, and ab-cd, respectively. Note that a, c and b, d do not necessarily share the same longitude/x-coordinate. See schematic:\n\n 0..............1 # fraction of distance Δab between a, b\n |< Δab >|\n\n 0^ a -------- o - b # anvil-shaped average of a, b, c, d at location x\n .Δy |\n . |\n .v x \n . |\n 1 c ------ o ---- d\n\n |< Δcd >|\n 0...............1 # fraction of distance Δcd between c, d\n\n^ fraction of distance Δy between a-b and c-d.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.average_on_poles-Union{Tuple{NF}, Tuple{AbstractGrid{NF}, Vector{<:UnitRange{<:Integer}}}} where NF<:Integer","page":"RingGrids","title":"SpeedyWeather.RingGrids.average_on_poles","text":"average_on_poles(\n A::AbstractGridArray{NF<:Integer, 1, Array{NF<:Integer, 1}},\n rings::Vector{<:UnitRange{<:Integer}}\n) -> Tuple{Any, Any}\n\n\nMethod for A::Abstract{T<:Integer} which rounds the averaged values to return the same number format NF.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.average_on_poles-Union{Tuple{NF}, Tuple{AbstractVector{NF}, Vector{<:UnitRange{<:Integer}}}} where NF<:AbstractFloat","page":"RingGrids","title":"SpeedyWeather.RingGrids.average_on_poles","text":"average_on_poles(\n A::AbstractArray{NF<:AbstractFloat, 1},\n rings::Vector{<:UnitRange{<:Integer}}\n) -> Tuple{Any, Any}\n\n\nComputes the average at the North and South pole from a given grid A and it's precomputed ring indices rings. The North pole average is an equally weighted average of all grid points on the northern-most ring. Similar for the South pole.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.check_inputs-NTuple{4, Any}","page":"RingGrids","title":"SpeedyWeather.RingGrids.check_inputs","text":"check_inputs(data, nlat_half, rings, Grid) -> Any\n\n\nTrue for data, nlat_half and rings that all match in size to construct a grid of type Grid.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.clenshaw_curtis_weights-Tuple{Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.clenshaw_curtis_weights","text":"clenshaw_curtis_weights(nlat_half::Integer) -> Any\n\n\nThe Clenshaw-Curtis weights for a Clenshaw grid (full or octahedral) of size nlathalf. Clenshaw-Curtis weights are of length nlat, i.e. a vector for every latitude ring, pole to pole. `sum(clenshawcurtisweights(nlathalf))is always2` as int0^π sin(x) dx = 2 (colatitudes), or equivalently int-pi/2^pi/2 cos(x) dx (latitudes).\n\nIntegration (and therefore the spectral transform) is exact (only rounding errors) when using Clenshaw grids provided that nlat >= 2(T + 1), meaning that a grid resolution of at least 128x64 (nlon x nlat) is sufficient for an exact transform with a T=31 spectral truncation.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.each_index_in_ring!-Tuple{Vector{<:UnitRange{<:Integer}}, Type{<:OctahedralGaussianArray}, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.each_index_in_ring!","text":"each_index_in_ring!(\n rings::Vector{<:UnitRange{<:Integer}},\n Grid::Type{<:OctahedralGaussianArray},\n nlat_half::Integer\n)\n\n\nprecompute a Vector{UnitRange{Int}} to index grid points on every ringj(elements of the vector) ofGridat resolutionnlat_half. Seeeachringandeachgrid` for efficient looping over grid points.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.each_index_in_ring-Tuple{Type{<:OctahedralGaussianArray}, Integer, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.each_index_in_ring","text":"each_index_in_ring(\n Grid::Type{<:OctahedralGaussianArray},\n j::Integer,\n nlat_half::Integer\n) -> Any\n\n\nUnitRange{Int} to index grid points on ring j of Grid at resolution nlat_half.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.each_index_in_ring-Tuple{Type{<:SpeedyWeather.RingGrids.AbstractFullGridArray}, Integer, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.each_index_in_ring","text":"each_index_in_ring(\n Grid::Type{<:SpeedyWeather.RingGrids.AbstractFullGridArray},\n j::Integer,\n nlat_half::Integer\n) -> Any\n\n\nUnitRange for every grid point of grid Grid of resolution nlat_half on ring j (j=1 is closest ring around north pole, j=nlat around south pole).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.each_index_in_ring-Union{Tuple{Grid}, Tuple{Grid, Integer}} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.each_index_in_ring","text":"each_index_in_ring(\n grid::AbstractGridArray,\n j::Integer\n) -> Any\n\n\nUnitRange to access data on grid grid on ring j.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.eachgrid-Tuple{AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.eachgrid","text":"eachgrid(grid::AbstractGridArray) -> Any\n\n\nCartesianIndices for the 2nd to last dimension of an AbstractGridArray, to be used like\n\nfor k in eachgrid(grid) for ring in eachring(grid) for ij in ring grid[ij, k]\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.eachgridpoint-Tuple{AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.eachgridpoint","text":"eachgridpoint(grid::AbstractGridArray) -> Base.OneTo\n\n\nUnitRange to access each horizontal grid point on grid grid. For a NxM (N horizontal grid points, M vertical layers) OneTo(N) is returned.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.eachgridpoint-Union{Tuple{Grid}, Tuple{Grid, Vararg{Grid}}} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.eachgridpoint","text":"eachgridpoint(\n grid1::AbstractGridArray,\n grids::Grid<:AbstractGridArray...\n) -> Base.OneTo\n\n\nLike eachgridpoint(::AbstractGridArray) but checks for equal size between input arguments first.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.eachring-Tuple{AbstractGridArray, Vararg{AbstractGridArray}}","page":"RingGrids","title":"SpeedyWeather.RingGrids.eachring","text":"eachring(\n grid1::AbstractGridArray,\n grids::AbstractGridArray...\n) -> Any\n\n\nSame as eachring(grid) but performs a bounds check to assess that all grids according to grids_match (non-parametric grid type, nlat_half and length).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.eachring-Tuple{AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.eachring","text":"eachring(grid::AbstractGridArray) -> Any\n\n\nVector{UnitRange} rings to loop over every ring of grid grid and then each grid point per ring. To be used like\n\nrings = eachring(grid)\nfor ring in rings\n for ij in ring\n grid[ij]\n\nAccesses precomputed grid.rings.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.eachring-Tuple{Type{<:AbstractGridArray}, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.eachring","text":"eachring(\n Grid::Type{<:AbstractGridArray},\n nlat_half::Integer\n) -> Any\n\n\nComputes the ring indices i0:i1 for start and end of every longitudinal point on a given ring j of Grid at resolution nlat_half. Used to loop over rings of a grid. These indices are also precomputed in every grid.rings.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.equal_area_weights-Tuple{Type{<:AbstractGridArray}, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.equal_area_weights","text":"equal_area_weights(\n Grid::Type{<:AbstractGridArray},\n nlat_half::Integer\n) -> Any\n\n\nThe equal-area weights used for the HEALPix grids (original or OctaHEALPix) of size nlathalf. The weights are of length nlat, i.e. a vector for every latitude ring, pole to pole. `sum(equalareaweights(nlathalf))is always2` as int0^π sin(x) dx = 2 (colatitudes), or equivalently int-pi/2^pi/2 cos(x) dx (latitudes). Integration (and therefore the spectral transform) is not exact with these grids but errors reduce for higher resolution.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.extrema_in-Tuple{AbstractVector, Real, Real}","page":"RingGrids","title":"SpeedyWeather.RingGrids.extrema_in","text":"extrema_in(v::AbstractVector, a::Real, b::Real) -> Any\n\n\nFor every element vᵢ in v does a<=vi<=b hold?\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.full_array_type-Tuple{AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.full_array_type","text":"full_array_type(grid::AbstractGridArray) -> Any\n\n\nFull grid array type for grid. Always returns the N-dimensional *Array not the two-dimensional (N=1) *Grid. For reduced grids the corresponding full grid that share the same latitudes.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.full_grid_type-Tuple{AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.full_grid_type","text":"full_grid_type(grid::AbstractGridArray) -> Any\n\n\nFull (horizontal) grid type for grid. Always returns the two-dimensional (N=1) *Grid type. For reduced grids the corresponding full grid that share the same latitudes.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.gaussian_weights-Tuple{Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.gaussian_weights","text":"gaussian_weights(nlat_half::Integer) -> Any\n\n\nThe Gaussian weights for a Gaussian grid (full or octahedral) of size nlathalf. Gaussian weights are of length nlat, i.e. a vector for every latitude ring, pole to pole. `sum(gaussianweights(nlathalf))is always2` as int0^π sin(x) dx = 2 (colatitudes), or equivalently int_-pi/2^pi/2 cos(x) dx (latitudes).\n\nIntegration (and therefore the spectral transform) is exact (only rounding errors) when using Gaussian grids provided that nlat >= 3(T + 1)/2, meaning that a grid resolution of at least 96x48 (nlon x nlat) is sufficient for an exact transform with a T=31 spectral truncation.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_colat-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_colat","text":"get_colat(grid::AbstractGridArray) -> Any\n\n\nColatitudes (radians) for meridional column (full grids only).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_colatlons-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_colatlons","text":"get_colatlons(grid::AbstractGridArray) -> Tuple{Any, Any}\n\n\nLatitudes (in radians, 0-π) and longitudes (0 - 2π) for every (horizontal) grid point in grid.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_lat-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_lat","text":"get_lat(grid::AbstractGridArray) -> Any\n\n\nLatitude (radians) for each ring in grid, north to south.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_latd-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_latd","text":"get_latd(grid::AbstractGridArray) -> Any\n\n\nLatitude (degrees) for each ring in grid, north to south.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_latdlonds-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_latdlonds","text":"get_latdlonds(grid::AbstractGridArray) -> Tuple{Any, Any}\n\n\nLatitudes (in degrees, -90˚-90˚N) and longitudes (0-360˚E) for every (horizontal) grid point in grid. Ordered 0-360˚E then north to south.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_latlons-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_latlons","text":"get_latlons(grid::AbstractGridArray) -> Tuple{Any, Any}\n\n\nLatitudes (in radians, 0-2π) and longitudes (-π/2 - π/2) for every (horizontal) grid point in grid.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_lon-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_lon","text":"get_lon(grid::AbstractGridArray) -> Any\n\n\nLongitude (radians) for meridional column (full grids only).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_lond-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_lond","text":"get_lond(grid::AbstractGridArray) -> Any\n\n\nLongitude (degrees) for meridional column (full grids only).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_nlat-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_nlat","text":"get_nlat(grid::AbstractGridArray) -> Any\n\n\nGet number of latitude rings, pole to pole.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_nlat_half-Tuple{AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_nlat_half","text":"get_nlat_half(grid::AbstractGridArray) -> Any\n\n\nResolution paraemeters nlat_half of a grid. Number of latitude rings on one hemisphere, Equator included.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_nlon_per_ring-Tuple{Type{<:HEALPixArray}, Integer, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_nlon_per_ring","text":"get_nlon_per_ring(\n Grid::Type{<:HEALPixArray},\n nlat_half::Integer,\n j::Integer\n) -> Any\n\n\nNumber of longitude points for ring j on Grid of resolution nlat_half.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_nlons-Tuple{Type{<:AbstractGridArray}, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_nlons","text":"get_nlons(\n Grid::Type{<:AbstractGridArray},\n nlat_half::Integer;\n both_hemispheres\n) -> Any\n\n\nReturns a vector nlons for the number of longitude points per latitude ring, north to south. Provide grid Grid and its resolution parameter nlat_half. For keyword argument both_hemispheres=false only the northern hemisphere (incl Equator) is returned.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_npoints-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_npoints","text":"get_npoints(grid::AbstractGridArray) -> Any\n\n\nTotal number of grid points in all dimensions of grid. Equivalent to length of the underlying data array.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.get_npoints2D-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.get_npoints2D","text":"get_npoints2D(grid::AbstractGridArray) -> Any\n\n\nNumber of grid points in the horizontal dimension only, even if grid is N-dimensional.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.grid_cell_average!-Tuple{AbstractGrid, SpeedyWeather.RingGrids.AbstractFullGrid}","page":"RingGrids","title":"SpeedyWeather.RingGrids.grid_cell_average!","text":"grid_cell_average!(\n output::AbstractGrid,\n input::SpeedyWeather.RingGrids.AbstractFullGrid\n) -> AbstractGrid\n\n\nAverages all grid points in input that are within one grid cell of output with coslat-weighting. The output grid cell boundaries are assumed to be rectangles spanning half way to adjacent longitude and latitude points.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.grid_cell_average-Tuple{Type{<:AbstractGrid}, Integer, SpeedyWeather.RingGrids.AbstractFullGrid}","page":"RingGrids","title":"SpeedyWeather.RingGrids.grid_cell_average","text":"grid_cell_average(\n Grid::Type{<:AbstractGrid},\n nlat_half::Integer,\n input::SpeedyWeather.RingGrids.AbstractFullGrid\n) -> Any\n\n\nAverages all grid points in input that are within one grid cell of output with coslat-weighting. The output grid cell boundaries are assumed to be rectangles spanning half way to adjacent longitude and latitude points.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.grids_match-Tuple{AbstractGridArray, AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.grids_match","text":"grids_match(\n A::AbstractGridArray,\n B::AbstractGridArray;\n horizontal_only,\n vertical_only\n) -> Any\n\n\nTrue if both A and B are of the same nonparametric grid type (e.g. OctahedralGaussianArray, regardless type parameter T or underyling array type ArrayType) and of same resolution (nlat_half) and total grid points (length). Sizes of (4,) and (4,1) would match for example, but (8,1) and (4,2) would not (nlat_half not identical).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.grids_match-Tuple{AbstractGridArray, Vararg{AbstractGridArray}}","page":"RingGrids","title":"SpeedyWeather.RingGrids.grids_match","text":"grids_match(\n A::AbstractGridArray,\n B::AbstractGridArray...;\n kwargs...\n) -> Any\n\n\nTrue if all grids A, B, C, ... provided as arguments match according to grids_match wrt to A (and therefore all).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.horizontal_grid_type-Tuple{AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.horizontal_grid_type","text":"horizontal_grid_type(grid::AbstractGridArray) -> Any\n\n\nThe two-dimensional (N=1) *Grid for grid, which can be an N-dimensional *GridArray.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.isdecreasing-Tuple{AbstractVector}","page":"RingGrids","title":"SpeedyWeather.RingGrids.isdecreasing","text":"isdecreasing(x::AbstractVector) -> Bool\n\n\nCheck whether elements of a vector v are strictly decreasing.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.isincreasing-Tuple{AbstractVector}","page":"RingGrids","title":"SpeedyWeather.RingGrids.isincreasing","text":"isincreasing(x::AbstractVector) -> Bool\n\n\nCheck whether elements of a vector v are strictly increasing.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.matrix_size-Tuple{Grid} where Grid<:AbstractGridArray","page":"RingGrids","title":"SpeedyWeather.RingGrids.matrix_size","text":"matrix_size(grid::AbstractGridArray) -> Tuple{Int64, Int64}\n\n\nSize of the matrix of the horizontal grid if representable as such (not all grids).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.nlat_odd-Tuple{AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.nlat_odd","text":"nlat_odd(grid::AbstractGridArray) -> Any\n\n\nTrue for a grid that has an odd number of latitude rings nlat (both hemispheres).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.nonparametric_type-Tuple{AbstractGridArray}","page":"RingGrids","title":"SpeedyWeather.RingGrids.nonparametric_type","text":"nonparametric_type(grid::AbstractGridArray) -> Any\n\n\nFor any instance of AbstractGridArray type its n-dimensional type (*Grid{T, N, ...} returns *Array) but without any parameters {T, N, ArrayType}\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.npoints_added_per_ring-Tuple{Type{<:OctahedralGaussianArray}}","page":"RingGrids","title":"SpeedyWeather.RingGrids.npoints_added_per_ring","text":"npoints_added_per_ring(\n _::Type{<:OctahedralGaussianArray}\n) -> Int64\n\n\n[EVEN MORE EXPERIMENTAL] number of longitude points added (removed) for every ring towards the Equator (on the southern hemisphere towards the south pole).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.npoints_pole-Tuple{Type{<:OctahedralGaussianArray}}","page":"RingGrids","title":"SpeedyWeather.RingGrids.npoints_pole","text":"npoints_pole(_::Type{<:OctahedralGaussianArray}) -> Int64\n\n\n[EXPERIMENTAL] additional number of longitude points on the first and last ring. Change to 0 to start with 4 points on the first ring.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.nside_healpix-Tuple{Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.nside_healpix","text":"nside_healpix(nlat_half::Integer) -> Any\n\n\nThe original Nside resolution parameter of the HEALPix grids. The number of grid points on one side of each (square) face. While we use nlat_half across all ring grids, this function translates this to Nside. Even nlat_half only.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.rotate_matrix_indices_180-Tuple{Integer, Integer, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.rotate_matrix_indices_180","text":"rotate_matrix_indices_180(\n i::Integer,\n j::Integer,\n s::Integer\n) -> Tuple{Any, Any}\n\n\nRotate indices i, j of a square matrix of size s x s by 180˚.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.rotate_matrix_indices_270-Tuple{Integer, Integer, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.rotate_matrix_indices_270","text":"rotate_matrix_indices_270(\n i::Integer,\n j::Integer,\n s::Integer\n) -> Tuple{Integer, Any}\n\n\nRotate indices i, j of a square matrix of size s x s anti-clockwise by 270˚.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.rotate_matrix_indices_90-Tuple{Integer, Integer, Integer}","page":"RingGrids","title":"SpeedyWeather.RingGrids.rotate_matrix_indices_90","text":"rotate_matrix_indices_90(\n i::Integer,\n j::Integer,\n s::Integer\n) -> Tuple{Any, Integer}\n\n\nRotate indices i, j of a square matrix of size s x s anti-clockwise by 90˚.\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.spherical_distance-Tuple{Type{<:SpeedyWeather.RingGrids.AbstractSphericalDistance}, Vararg{Any}}","page":"RingGrids","title":"SpeedyWeather.RingGrids.spherical_distance","text":"spherical_distance(\n Formula::Type{<:SpeedyWeather.RingGrids.AbstractSphericalDistance},\n args...;\n kwargs...\n) -> Any\n\n\nSpherical distance, or great-circle distance, between two points lonlat1 and lonlat2 using the Formula (default Haversine).\n\n\n\n\n\n","category":"method"},{"location":"ringgrids/#SpeedyWeather.RingGrids.whichring-Tuple{Integer, Vector{UnitRange{Int64}}}","page":"RingGrids","title":"SpeedyWeather.RingGrids.whichring","text":"whichring(\n ij::Integer,\n rings::Vector{UnitRange{Int64}}\n) -> Int64\n\n\nObtain ring index j from gridpoint ij and rings describing rind indices as obtained from eachring(::Grid)\n\n\n\n\n\n","category":"method"},{"location":"structure/#Tree-structure","page":"Tree structure","title":"Tree structure","text":"","category":"section"},{"location":"structure/","page":"Tree structure","title":"Tree structure","text":"At the top of SpeedyWeather's type tree sits the Simulation, containing variables and model, which in itself contains model components with their own fields and so on. (Note that we are talking about the structure of structs within structs not the type hierarchy as defined by subtyping abstract types.) This can quickly get complicated with a lot of nested structs. The following is to give users a better overview of how simulation, variables and model are structured within SpeedyWeather. Many types in SpeedyWeather have extended Julia's show function to give you an overview of its contents, e.g. a clock::Clock is printed as","category":"page"},{"location":"structure/","page":"Tree structure","title":"Tree structure","text":"using SpeedyWeather\nclock = Clock()","category":"page"},{"location":"structure/","page":"Tree structure","title":"Tree structure","text":"illustrating the fields within a clock, their types and (unless they are an array) also their values. For structs within structs however, this information, would not be printed by default. You could use Julia's autocomplete like clock.
