\section{Advection} \label{sec:adv}
Advection is a fluid flow transporting something with it as it flows. This can be temperature, gas, solids or other fluids. In our case we will be looking at temperature.

\subsection{Thermal Diffusion}
As of this time, what you notice if you run the model is that the winds only get stronger and stronger (and the model is hence blowing up, which means that the numbers increase so dramatically 
that it is no longer realistic). This is because there is no link yet between the velocities of the atmosphere and the temperature. Currently, any air movement does not affect the temperature 
in the atmosphere of our model while it does in reality. So we need to change some calculations to account for that. Thermal diffusion helps with spreading out the temperatures and tempering 
the winds a bit.

The diffusion equation, as written in \autoref{eq:diffusion}, describes how the temperature spreads out over time\cite{diffusion}. The symbols in the equation represent:

\begin{itemize}
    \item $u$: A vector consisting out of 4 elements: $x, y, z, t$. $x, y, z$ are the local coordinates and $t$ is time.
    \item $\alpha$: The thermal diffusivity constant.
    \item $\nabla^2$: The Laplace operator, more information in \autoref{sec:laplace}.
    \item $\bar{u}$: The time derivative of $u$, or in symbols $\frac{\delta u}{\delta t}$.
\end{itemize}

\begin{equation}
    \bar{u} = \alpha \nabla^2 u
    \label{eq:diffusion}
\end{equation}

Now to get this into code we need the following algorithms \autoref{alg:laplacian} and \autoref{alg:diffusion}. \autoref{alg:laplacian} implements the laplacian operator, whereas 
\autoref{alg:diffusion} implements the diffusion calculations.  $\nabla^2$ in \autoref{alg:diffusion} represents the call to \autoref{alg:laplacian}.

\begin{algorithm}
    $T_a \leftarrow T_a + \delta t \alpha_a \nabla^2(T_a)$ \;
    $T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
    \caption{The main calculations for calculating the effects of diffusion}
    \label{alg:diffusion}
\end{algorithm}

\subsection{Adding in Advection}
With thermal diffusion in place, the temperature will spread out a bit, however air is not transported yet. This means that the winds we simulate are not actually moving any air. Advection is
going to change that. The advection equation is shown in \autoref{eq:advection}. The symbols are:

\begin{itemize}
    \item $\psi$: What is carried along (in our case temperature, $K$).
    \item $t$: The time ($s$).
    \item $u$: The fluid velocity vector ($ms^{-1}$).
    \item $\nabla$: The divergence operator (as explained in \autoref{sec:laplace}).
\end{itemize}

\begin{equation}
    \frac{\delta \psi}{\delta t} + \nabla \cdot (\psi u) = 0
    \label{eq:advection}
\end{equation}

With the divergence functon defined in \autoref{alg:divergence}, we now need to adjust \autoref{alg:diffusion} to incorporate this effect. The resulting algorithm can be found in 
\autoref{alg:advection}. Here $\nabla$ represents the function call to \autoref{alg:divergence}.

\begin{algorithm}
    $T_{add} \leftarrow T_a + \delta t \alpha_a \nabla^2(T_a) + \nabla(T_a)$ \;
    $T_a \leftarrow T_a + T_{add}[5:-5, :] \text{ //Only add } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + 5, nlat - 5]$. \;
    $T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
    \caption{The main calculations for calculating the effects of diffusion}
    \label{alg:advection}
\end{algorithm}

Now that we have the air moving, we also need to account for the moving of the density. This is because moving air to a certain place will change the air density at that place if the air at that 
place does not move away at the same rate. Say we are moving air to $x$ at $y \ ms^{-1}$. If air at $x$ moves at a rate $z \ ms^{-1}$ and $z \neq y$ then the air density at $x$ will change.
The equation we will need for that is the mass continuity equation as shown in \autoref{eq:mass continuity} \cite{masscontinue}.

\begin{equation}
    \frac{\delta \rho}{\delta t} + \nabla \cdot (\rho v) = 0
    \label{eq:mass continuity}
\end{equation}

Using this equation means that we will no longer assume that the atmosphere is incompressible. Therefore we need to change a few things in the code. First we need to change the $\rho$ in 
\autoref{alg:stream3}. Since $\rho$ is no longer constant we need to access the right value of $\rho$ by specifying the indices. So $\rho$ will change to $\rho[lat, lon]$. Furthermore we need
to calculate $\rho$ after the movement of air has taken place, so we need to change \autoref{alg:advection} as well to include the calculations for $\rho$. The new version can be found in 
\autoref{alg:advectionv2}. Again the $\nabla$ represents the call to \autoref{alg:divergence}.


\begin{algorithm}
    $T_{add} \leftarrow T_a + \delta t \alpha_a \nabla^2(T_a) + \nabla(T_a)$ \;
    $T_a \leftarrow T_a + T_{add}[5:-5, :] \text{ //Only add } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + 5, nlat - 5]$. \;
    $\rho \leftarrow \rho + \delta t \nabla \rho$ \;
    $T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
    \caption{The main calculations for calculating the effects of diffusion}
    \label{alg:advectionv2}
\end{algorithm}

Currently the advection does not work like it should. This is probably due to boundary issues, where we get too close to the poles and it starts freaking out there \cite{simon}. So to fix this 
we are going to define boundaries and assume that the advection only works within those boundaries. We only let it change by half of the values. The changes are incorporated in 
\autoref{alg:advectionfix}. The reason why we mention this seperately, in contrast to the other fixes that we have incorporated throughout the manual already, is the accompanying change with the 
boundary. 

\begin{algorithm}
        $T_{add} \leftarrow T_a + \delta t \alpha_a \nabla^2(T_a) + \nabla(T_a)$ \;
        $T_a \leftarrow T_a - 0.5T_{add}[adv\_bound:-adv\_boun, :] \text{ //Only subtract } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$. \;
        $\rho[adv\_boun: -adv\_boun, :] \leftarrow \rho - 0.5(\delta t \nabla \rho) \text{ //Only change the density for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$ \;
        $T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
    \caption{The main calculations for calculating the effects of diffusion}
    \label{alg:advectionfix}
\end{algorithm}

\subsection{Layers, layers and layers}
With the atmospheric layers, and all matrices that have an extra dimension to account for it, we need to add the correct indices to the advection algorithm \autoref{alg:advectionfix}. Let us 
add it, with \autoref{alg:advection layer} as a result. Here the ':' means all indices of the 3 dimensional matrix. 

\begin{algorithm}
    $T_{add} \leftarrow T_a + \delta t \alpha_a \nabla^2(T_a) + \nabla(T_a)$ \;
    $T_a \leftarrow T_a - 0.5T_{add}[adv\_boun:-adv\_boun, :, :] \text{ //Only subtract } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$. \;
    $\rho[adv\_boun: -adv\_boun, :, :] \leftarrow \rho - 0.5(\delta t \nabla \rho) \text{ //Only change the density for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$ \;
    $T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
    \caption{The main calculations for calculating the effects of diffusion}
    \label{alg:advection layer}
\end{algorithm}

First thing to mention is that vertical advection is still broken. Why? Because the gradient in the $z$ direction is broken. This is due to finite differencing on an exponential function. The way
we calculate the difference from one layer to the other is by differencing them (subtracting) which is always finite. Therefore we always get some inaccuracies. Usually that is fine, but with an 
exponential function the differences, you guessed it, become exponentially wrong. As such, the function would eventually be so far off that the model would blow up. So we need to fix it. To 
prevent a blow up, we have disabled the call to the gradient $z$ function in \autoref{alg:divergence}. This ensures that the horizontal bits still work, but the vertical stuff does not.
As always, we will try to fix this in a future stream. 