mirror of https://github.com/Askill/claude.git
231 lines
15 KiB
TeX
231 lines
15 KiB
TeX
\section{Let's Get the Atmosphere Moving}
|
|
In its current state, CLaUDE has a static planet. This means that the planet remains in place and does not move. However we know that planets move in orbit and more importantly, spin around
|
|
themselves. But before we start adding layers, let's talk about a term you will hear more often: numerical instability.
|
|
|
|
Numerical instability occurs when you first run the model. This is due to the nature of the equations. Nearly all equations are continuous, which means that they are always at work. However
|
|
when you start the model, the equations were not at work yet. It is as if you suddenly give a random meteor an atmosphere, place it in orbit around a star and don't touch i for a bit. You will
|
|
see that the whole system oscilates wildly as it adjusts to the sudden changes and eventually it will stabilise. Another term you might encounter is blow up, this occurs when when the model
|
|
suddenly no longer behaves like it should. This is most likely caused by mistakes in the code or incorrect paramter initialisation. Be wary of the existence of both factors, and do not dismiss
|
|
a model if it behaves weirdly as it has just started up.
|
|
|
|
\subsection{Equation of State and the Incompressible Atmosphere}
|
|
The equation of state relates one or more variables in a dynamical system (like the atmosphee) to another. The most common equation of state in the atmosphere is the ideal gas equation as
|
|
described by \autoref{eq:ideal gas} \cite{idealGas}. The symbols in that equation represent:
|
|
|
|
\begin{itemize}
|
|
\item $p$: The gas pressure ($Pa$).
|
|
\item $V$: The volume of the gas ($m^3$).
|
|
\item $n$: The amount of moles\footnote{Mole is the amount of particles ($6.02214076 \cdot 10^{23}$) in a substance, where the average weight of one mole of particles in grams is about the
|
|
same as the weight of one particle in atomic mass units ($u$)\cite{mole}} in the gas.
|
|
\item $R$: The Gas constant, $8.3144621$ ($J(mol)^{-1}K$) \cite{idealGas}.
|
|
\item $T$: The temperature opf the gas ($K$).
|
|
\end{itemize}
|
|
|
|
If we divide everything in \autoref{eq:ideal gas} by $V$ and set it to be unit (in this case, set it to be exactly $1 m^3$) we can add in the molar mass in both the top and bottom parts of the
|
|
division as show in \autoref{eq:gas unit}. We can then replace $\frac{nm}{V}$ by $\rho$ the density of the gas ($kgm^{-3}$) and $\frac{R}{m}$ by $R_s$ the specific gas constant (gas constant that varies per
|
|
gas in $J(mol)^{-1}K$) as shown in \autoref{eq:state gas}. the resulting equation is the equation of state that you get that most atmospheric physicists use when talking about the atmosphere \cite{simon}.
|
|
|
|
\begin{subequations}
|
|
\begin{equation}
|
|
pV = nRT
|
|
\label{eq:ideal gas}
|
|
\end{equation}
|
|
\begin{equation}
|
|
p = \frac{nR}{V}T = \frac{nmR}{Vm}T
|
|
\label{eq:gas unit}
|
|
\end{equation}
|
|
\begin{equation}
|
|
p = \rho R_sT
|
|
\label{eq:state gas}
|
|
\end{equation}
|
|
\end{subequations}
|
|
|
|
The pressure is quite important, as air moves from a high pressure point to a low pressure point. So if we know the density and the temperature, then we know the pressure and we can work out
|
|
where the air will be moving to (i.e. how the wind will blow). In our current model, we know the atmospheric temperature but we do not know the density. For simplicities sake, we will now assume
|
|
that the atmosphere is Incompressible, meaning that we have a constant density. Obviously we know that air can be compressed and hence our atmosphere can be compressed too but that is not
|
|
important enough to account for yet, especially considering the current complexity of our model.
|
|
|
|
The code that corresponds to this is quite simple, the only change that we need to make in \autoref{eq:state gas} is that we need to replace $T$ by $T_a$, the temperature of the atmosphere. As
|
|
$T_a$ is a matrix (known to programmers as a double array), $p$ will be a matrix as well. Now we only need to fill in some values. $\rho = 1.2$\cite{densityAir}, $R_s = 287$\cite{specificGasConstantAir}.
|
|
|
|
\subsection{The Primitive Equations and Geostrophy}
|
|
The primitive equations (also known as the momentum equations) is what makes the air move. It is actually kind of an injoke between physicists as they are called the primitive equations but
|
|
actually look quite complicated (and it says $fu$ at the end! \cite{simon}). The primitive equations are a set of equations dictating the direction in the $u$ and $v$ directions as shown in
|
|
\autoref{eq:primitive u} and \autoref{eq:primitive v}. We can make the equations simpler by using and approximation called geostrophy which means that we have no vertical motion, such that the
|
|
terms with $\omega$ in \autoref{eq:primitive u} and \autoref{eq:primitive v} become 0. We also assume that we are in a steady state, i.e. there is no acceleration which in turn means that the
|
|
whole middle part of the equations are $0$. Hence we are left with \autoref{eq:primitive u final} and \autoref{eq:primitive v final}.
|
|
|
|
\begin{subequations}
|
|
\begin{equation}
|
|
\frac{du}{dt} = \frac{\delta u}{\delta t} + u\frac{\delta u}{ \delta x} + v\frac{\delta u}{\delta v} + \omega\frac{\delta u}{\delta p} = -\frac{\delta \Phi}{\delta x} + fv
|
|
\label{eq:primitive u}
|
|
\end{equation}
|
|
\begin{equation}
|
|
\frac{dv}{dt} = \frac{\delta v}{\delta t} + u\frac{\delta v}{ \delta x} + v\frac{\delta v}{\delta v} + \omega\frac{\delta v}{\delta p} = -\frac{\delta \Phi}{\delta y} - fu
|
|
\label{eq:primitive v}
|
|
\end{equation}
|
|
|
|
\begin{equation}
|
|
0 = -\frac{\delta \Phi}{\delta x} + fv
|
|
\label{eq:primitive u final}
|
|
\end{equation}
|
|
\begin{equation}
|
|
0 = -\frac{\delta \Phi}{\delta y} - fu
|
|
\label{eq:primitive v final}
|
|
\end{equation}
|
|
\end{subequations}
|
|
|
|
\autoref{eq:primitive u final} can be split up into to parts, the $\frac{\delta \Phi}{\delta x}$ part (the gradient force) and the $fv$ part (the coriolis force). The same applies to
|
|
\autoref{eq:primitive v final}. Effectively we have a balance between the gradient and the coriolis force as shown in \autoref{eq:pu simple} and \autoref{eq:pv simple}. The symbols in both of
|
|
these equations are:
|
|
|
|
\begin{itemize}
|
|
\item $\Phi$: The geopotential, potential (more explanation in \autoref{sec:potential}) of the planet's gravity field ($Jkg^{-1}$).
|
|
\item $x$: The change in the East direction along the planet surface ($m$).
|
|
\item $y$: The change in the North direction along the planet surface ($m$).
|
|
\item $f$: The coriolis parameter as described by \autoref{eq:coriolis}, where $\Omega$ is the rotation rate of the planet (for Earth $7.2921 \cdot 10^{-5}$) ($rad \ s^{-1}$) and $\theta$ is the
|
|
latitude \cite{coriolis}.
|
|
\item $u$: The velocity in the latitude ($ms^{-1}$).
|
|
\item $v$: The velocity in the longitude ($ms^{-1}$).
|
|
\end{itemize}
|
|
|
|
\begin{subequations}
|
|
\begin{equation}
|
|
f = 2\Omega\sin(\theta)
|
|
\label{eq:coriolis}
|
|
\end{equation}
|
|
\begin{equation}
|
|
\frac{\delta \Phi}{\delta x} = fv
|
|
\label{eq:pu simple}
|
|
\end{equation}
|
|
\begin{equation}
|
|
\frac{\delta \Phi}{\delta y} = -fu
|
|
\label{eq:pv simple}
|
|
\end{equation}
|
|
\begin{equation}
|
|
\frac{\delta p}{\rho \delta x} = fv
|
|
\label{eq:pu simple final}
|
|
\end{equation}
|
|
\begin{equation}
|
|
\frac{\delta p}{\rho \delta y} = -fu
|
|
\label{eq:pv simple final}
|
|
\end{equation}
|
|
\end{subequations}
|
|
|
|
Since we want to know how the atmosphere moves, we want to get the v and u components of the velocity vector (since $v$ and $u$ are the veolicites in longitude and latitude, if we combine them in a
|
|
vector we get the direction of the overall velocity). So it is time to start coding and calculating! If we look back at \autoref{alg:stream1v2}, we can see that we already have a double for loop.
|
|
In computer science, having multiple loops is generally considered a bad coding practice as you usually can just reuse the indices of the already existing loop, so you do not need to create a new
|
|
one. However this is a special case, since we are calculating new temperatures in the double for loop. If we then also would start to calculate the velocities then we would use new information
|
|
and old information at the same time. Since at index $i - 1$ the new temperature has already been calculated, but at the index $i + 1$ the old one is still there. So in order to fix that we need
|
|
a second double for loop to ensure that we always use the new temperatures. We display this specific loop in \autoref{alg:stream2}. Do note that everything in \autoref{alg:stream1v2} is still
|
|
defined and can still be used, but since we want to focus on the new code, we leave out the old code to keep it concise and to prevent clutter.
|
|
|
|
\begin{algorithm}[hbt]
|
|
\SetAlgoLined
|
|
|
|
\While{\texttt{TRUE}}{
|
|
\For{$lat \in [-nlat, nlat]$}{
|
|
\For{$lon \in [0, nlon]$}{
|
|
$u[lat, lon] \leftarrow -\frac{p[lat + 1, lon] - p[lat - 1, lon]}{\delta y} \cdot \frac{1}{f[lat]\rho}$ \;
|
|
$v[lat, lon] \leftarrow \frac{p[lat, lon + 1] - p[lat, lon - 1]}{\delta x[lat]} \cdot \frac{1}{f[lat]\rho}$ \;
|
|
}
|
|
}
|
|
}
|
|
\caption{The main loop of the velocity of the atmosphere calculations}
|
|
\label{alg:stream2}
|
|
\end{algorithm}
|
|
|
|
The gradient calculation is done in \autoref{alg:gradient}. For this to work, we need the circumference of the planet. Herefore we need to assume that the planet is a sphere. While that is not
|
|
technically true, it makes little difference in practice and is good enough for our model. The equation for the circumference can be found in \autoref{eq:circumference} \cite{circumference},
|
|
where $r$ is the radius of the planet. Here we also use the f-plane approximation, where the coriolis paramter has one value for the northern hemisphere and one value for the southern hemisphere
|
|
\cite{fplane}.
|
|
|
|
\begin{equation}
|
|
2 \pi r
|
|
\label{eq:circumference}
|
|
\end{equation}
|
|
|
|
\begin{algorithm}
|
|
\SetAlgoLined
|
|
$C \leftarrow 2\pi R$ \;
|
|
$\delta y \leftarrow \frac{C}{nlat}$ \;
|
|
|
|
\For{$lat \in [-nlat, nlat]$}{
|
|
$\delta x[lat] \leftarrow \delta y \cos(lat \cdot \frac{\pi}{180})$ \;
|
|
|
|
\eIf{$lat < 0$}{
|
|
$f[lat] \leftarrow -10^{-4}$ \;
|
|
}{
|
|
$f[lat] \leftarrow 10^{-4}$ \;
|
|
}
|
|
}
|
|
\caption{Calculating the gradient $\delta x$}
|
|
\label{alg:gradient}
|
|
\end{algorithm}
|
|
|
|
Because of the geometry of the planet and the construction of the longitude latitude grid, we run into some problems when calculating the gradient. Since the planet is not flat ("controversial
|
|
I know"\cite{simon}) whenever we reach the end of the longitude we need to loop around to get to the right spot to calculate the gradients (as the planet does not stop at the end of the
|
|
longitude line but loops around). So to fix that we use the modulus (mod) function which does the looping for us if we exceed the grid's boundaries. We do haveanother problem though, the poles.
|
|
As the latitude grows closer to the poles, they are converging on the center point of the pole. Looping around there is much more difficult so to fix it, we just do not consider that center
|
|
point in the main loop. The changed algorithm can be found in \autoref{alg:stream2v2}
|
|
|
|
\begin{algorithm}[hbt]
|
|
\SetAlgoLined
|
|
|
|
\While{\texttt{TRUE}}{
|
|
\For{$lat \in [-nlat + 1, nlat - 1]$}{
|
|
\For{$lon \in [0, nlon]$}{
|
|
$u[lat, lon] \leftarrow -\frac{p[(lat + 1) \text{ mod } nlat, lon] - p[(lat -1) \text{ mod } nlat, lon]}{\delta y} \cdot \frac{1}{f[lat]\rho}$ \;
|
|
$v[lat, lon] \leftarrow \frac{p[lat, (lon + 1) \text{ mod } nlon] - p[lat, (lon -1) \text{ mod } nlon]}{\delta x[lat]} \cdot \frac{1}{f[lat]\rho}$ \;
|
|
}
|
|
}
|
|
}
|
|
\caption{The main loop of the velocity of the atmosphere calculations}
|
|
\label{alg:stream2v2}
|
|
\end{algorithm}
|
|
|
|
Do note that the pressure calculation is done between the temperature calculation in \autoref{alg:stream1v2} and the $u, v$ calculations in \autoref{alg:stream2v2}. At this point our model shows
|
|
a symmetric vortex around the sun that moves with the sun. This is not very realistic as you usually have convection and air flowing from warm to cold, but we do not have that complexity yet
|
|
(due to our single layer atmosphere).
|
|
|
|
\subsection{Introducing an Ocean}
|
|
Now we want to introduce an ocean, because most of the Earth is covered by oceans it plays quite an important role in atmospheric physics. To do this we need a new concept called albedo. Albedo
|
|
is basically the reflectiveness of a material (in our case the planet's surface) \cite{albedo}. The average albedo of the Earth is about 0.3. Now to add an ocean to the grid, we define a few
|
|
areas where the albedo differs. Where you do this does not really matter for the current complexity. Defining the oceans is as easy as hardcoding (what we computer scientists refer to when
|
|
setting parts of an array to be a specific value, where if you want to change the value you need to change it everywhere instead of doing it in a variable) the albedo value for the specific
|
|
regions as we do in \autoref{alg:albedo}. Water also takes longer to warm up, so let us change the specific heat capacity ($C_p$ in \autoref{alg:stream1v2}) from a constant to an array. The new
|
|
$C_p$ can also be found in \autoref{alg:albedo}, where we have made the specific heat capacity of water one order of magnitude (i.e. $10$ times) larger.
|
|
|
|
\begin{algorithm}[hbt]
|
|
$a \leftarrow 0.5$ \;
|
|
$a[5-55, 9-20] \leftarrow 0.2$ \;
|
|
$a[23-50, 45-70] \leftarrow 0.2$ \;
|
|
$a[2-30, 85-110] \leftarrow 0.2$ \;
|
|
|
|
$C_p \leftarrow 10^7$ \;
|
|
$C_p[5-55, 9-20] \leftarrow 10^8$ \;
|
|
$C_p[23-50, 45-70] \leftarrow 10^8$ \;
|
|
$C_p[2-30, 85-110] \leftarrow 10^8$ \;
|
|
\caption{Defining the oceans}
|
|
\label{alg:albedo}
|
|
\end{algorithm}
|
|
|
|
Now that we have that defined, we need to adjust the main loop of the program (\autoref{alg:stream1v2}). For clarity, all the defined constants have been left out. We need to add albedo into the
|
|
equation and change $C_p$ from a constant to an array. The algorithm after these changes can be found in \autoref{alg:stream2v3}. We multiply by $1 - a$ since albedo represents how much energy is
|
|
reflected instead of absorbed, where we need the amount that is absorbed which is exactly equal to $1$ minus the amount that is reflected.
|
|
|
|
\begin{algorithm}[hbt]
|
|
\SetAlgoLined
|
|
|
|
\While{\texttt{TRUE}}{
|
|
\For{$lat \in [-nlat, nlat]$}{
|
|
\For{$lon \in [0, nlot]$}{
|
|
$T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t ((1 - a[lat, lon])S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p[lat, lon]}$ \;
|
|
$T_a[lat, lon] \leftarrow T_a[lat, lon] + \frac{\delta t (\sigma (T_p[lat, lon])^4 - 2\epsilon\sigma (T_a[lat, lon])^4)}{C_a}$ \;
|
|
$t \leftarrow t + \delta t$ \;
|
|
}
|
|
}
|
|
}
|
|
\caption{The main loop of the temperature calculations}
|
|
\label{alg:stream2v3}
|
|
\end{algorithm} |