askill
/
claude
mirror of https://github.com/Askill/claude.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Updated the manual to incorporate the latest stream and putting the documentation source files on the repo

This commit is contained in:
TechWizzart 2020-07-16 18:40:06 +02:00
parent 249c6e1d4e
commit 147d42d418
No known key found for this signature in database
GPG Key ID: 151E66B6CEF55F2E
12 changed files with 1609 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
CLAuDE NOM.pdf
View File

Binary file not shown.

101
Documentation Source/CLAuDE.tex Normal file
View File

@ -0,0 +1,101 @@
\documentclass{article}
\usepackage{hyperref}
\usepackage{cite}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage[ruled]{algorithm2e}
\usepackage[margin=1in]{geometry}
\usepackage{appendix}
\usepackage{float}
\usepackage{verbatim}
\setcounter{section}{-1}
\title{CLimate Analysis using Digital Estimations Non-Offical Manual (CLAuDE NOM)}
\author{TechWizard}
\date{\today}
\begin{document}
\maketitle
\tableofcontents
\section{Introduction}
The CLimate Analysis using Digital Estimations model is a simplified planetary climate model. It will be used to educate people on how climate physics works and to experiment with different
parameters and see how much influence a tiny change can have (like for instance the rotation rate of the planet around its axis). It is built to be accessible to and runnable by everyone,
whether they have a super computer or a dated laptop. The model is written in Python and written during the weekly streams of Dr. Simon Clark \cite{twitch}. Each subsequent section starts with a
number, this number indicates which coding stream corresponds to that section. This does not only make it easier to cross reference but if the explanation is unclear or you just want to watch
the stream about that specific topic, you know which stream to watch. There is a useful playlist on Simon's Twitch which has all the streams without ad breaks or interruptions \cite{playlist}.
This manual will provide an overview of the formulae used and will explain aspects of these formulae. For each equation each symbol will be explained what it is. In such an explanation, the units
will be presented in SI units \cite{SI} between brackets like: $T$: The temperature of the planet ($K$). Which indicates that $T$ is the temperature of the planet in degrees Kelvin. If you need
to relate SI units to your preferred system of units, please refer to the internet for help with that. There are great calculators online where you only need to plug in a number and select the
right units.
Within this manual we will not concern ourselves with plotting the data, instead we focus on the physics side of things and translating formular into code. If you are interested in how the
plotting of the data works, or how loading and saving data works, please refer to the relevant stream on Simon's Twitch page \cite{twitch}.
This manual is for the toy model, which is as of now still in development. Therefore this manual will be in chronological order, explaining everything in the same order as it has been done.
There are plans to eventually modularise the whole model into seperate parts that can be extended by the community. When that will hit development, a new manual for that version of the model
will be made that treats things per topic instead of chronological order.
One important thing to note, the layout may change significantly when new sections are added. This is due to the amount of code that is added/changed. If a lot of code changes, a lot of so called
'algorithm' blocks are present which have different placement rules than just plain text. Therefore it may occur that an algorithm is referenced even though it is one or two pages later. This is
a pain to fix and if something later on changes the whole layout may be messed up again and is a pain to fix again. Hence I opt to let \LaTeX (the software/typeset language used to create this
manual) figure out the placement of the algorithm blocks, which may or may not be in the right places.
\input{Streams/Stream1.tex}
\input{Streams/Stream2.tex}
\input{Streams/Stream3.tex}
\input{Streams/Stream4.tex}
\input{Streams/Stream5.tex}
\input{Streams/Stream6.tex}
\input{Streams/Stream7.tex}
\newpage
\input{Streams/TTNMETAF.tex}
\newpage
\bibliography{references}
\bibliographystyle{plain}
\begin{comment}
\subsection{Hydrostatic balance}
Hydrostatic balance is the relation between the pressure at the surface of the planet and the pressure higher up in the atmosphere \cite{hydrostatic}. This all works under the assumption that
the vertical acceleration is relatively small in magnitude compared to the gravity of the planet, which it almost always is when looking at the atmosphere or the ocean. The corresponding
formula is described in \autoref{eq:hydrostatic}. The symbols here represent:
\begin{itemize}
\item $\delta p$: The change in the atmospheric pressure ($Pa$).
\item $\delta z$: The change in altitude or height ($m$).
\item $\rho$: The density of the atmosphere ($kgm^{-3}$).
\item $g$: The gravitational acceleration, on Earth this is $9.81$ ($ms^{-2}$).
\end{itemize}
\begin{equation}
\frac{\delta p}{\delta z} = -\rho g
\label{eq:hydrostatic}
\end{equation}
Now we need to translate this into code, luckily this is quite easy as shown in \autoref{alg:hydrostatic}. We first calculate the difference between the current mean pressure and the mean
pressure we get from \autoref{eq:hydrostatic}. Then we add the difference to the pressure, effectively correcting the average pressure for the current atmospheric layer to what it is supposed to
be. By correcting it with the mean difference we still have the local variation in pressure while on average we are closer to the true value we calculate in \autoref{eq:hydrostatic}.
\begin{algorithm}
\For{$alt \in [1, nlevels]$}{
$correction \leftarrow mean(p[:, :, alt - 1]) - mean(\rho[:, :, alt])(heights[alt] - heights[alt - 1])g$ \;
$p[:, :, alt] \leftarrow \rho[:, :, alt] + correction$ \;
}
\caption{Calculating the effects of Hydrostatic Balance on the average pressure}
\label{alg:hydrostatic}
\end{algorithm}
\end{comment}
\end{document}

29
Documentation Source/Streams/Stream.aux Normal file
View File

@ -0,0 +1,29 @@
\relax
\providecommand\hyper@newdestlabel[2]{}
\@setckpt{Streams/"Stream}{
\setcounter{page}{2}
\setcounter{equation}{0}
\setcounter{enumi}{0}
\setcounter{enumii}{0}
\setcounter{enumiii}{0}
\setcounter{enumiv}{0}
\setcounter{footnote}{0}
\setcounter{mpfootnote}{0}
\setcounter{part}{0}
\setcounter{section}{0}
\setcounter{subsection}{0}
\setcounter{subsubsection}{0}
\setcounter{paragraph}{0}
\setcounter{subparagraph}{0}
\setcounter{figure}{0}
\setcounter{table}{0}
\setcounter{Item}{0}
\setcounter{Hfootnote}{0}
\setcounter{bookmark@seq@number}{1}
\setcounter{parentequation}{0}
\setcounter{AlgoLine}{0}
\setcounter{algocfline}{0}
\setcounter{algocfproc}{0}
\setcounter{algocf}{0}
\setcounter{section@level}{1}
}

228
Documentation Source/Streams/Stream1.tex Normal file
View File

@ -0,0 +1,228 @@
\section{The Beginning}
\subsection{The First Law of Thermodynamics and the Stefan-Boltzmann Equation}
The beginning of CLAuDE is based upon one of the most important laws of physics: "Energy is neither created nor destroyed, only changed from one form to another." In otherwords, if energy goes into an object it must
equal the outflowing energy plus the change of internal energy. This is captured in Stefan-Boltzmann's law (\autoref{eq:stefan-boltzmann}) \cite{stefan-boltzmann}.
Here we assume that the planet is a black body, i.e. it absorbs all radiation (energy waves, some waves are visible like light, others are invisible like radio signals) on all wavelengths.
In \autoref{eq:stefan-boltzmann} the symbols are:
\begin{itemize}
\item $S$: The energy that reaches the top of the atmosphere, coming from the sun or a similar star, per meters squared $Jm^{-2}$. This is also called the insolation.
\item $\sigma$: The Stefan-Boltzmann constant, $5.670373 \cdot 10^-8 \ (Wm^{-2}K^{-4})$ \cite{stefan-boltzmann}.
\item $T$: The temperature of the planet ($K$).
\end{itemize}
Technically speaking \autoref{eq:stefan-boltzmann} is incorrect, as there is a mismatch in units. However, that is corrected in \autoref{eq:basis sphere final} so there is no need to worry about
it just yet. The energy difference between the energu that reaches the atmosphere and the temperature of the planet must be equal to the change in temperature of the planet, which is written in
\autoref{eq:sb rewritten}. The symbols on the right hand side are:
\begin{itemize}
\item $\Delta U$: The change of internal energy ($J$) \cite{thermo1}.
\item $C$: The specific heat capacity of the object, i.e. how much energy is required to heat the object by one degree Kelvin ($\frac{J}{K}$).
\item $\Delta T$: The change in temperature ($K$).
\end{itemize}
We want to know the change of temperature $\Delta T$, so we rewrite the equation into \autoref{eq:sb rewritten2}. Here we added the $\delta t$ term to account for the time difference (or time step). This is needed as
we need an interval to calculate the difference in temperature over. Also we needed to make the units match, and by adding this time step the units all match up perfectly.
\begin{subequations}
\begin{equation}
S = SB = \sigma T^4
\label{eq:stefan-boltzmann}
\end{equation}
\begin{equation}
S - \sigma T^4 = \Delta U = C \Delta T
\label{eq:sb rewritten}
\end{equation}
\begin{equation}
\Delta T = \frac{\delta t(S - \sigma T^4)}{C}
\label{eq:sb rewritten2}
\end{equation}
\label{eq:basis}
\end{subequations}
The set of equations in \autoref{eq:basis} form the basis of the temperature exchange of the planet. However two crucial aspects are missing. Only half of the planet will be receiving light from
the sun at once, and the planet is a sphere. So we need to account for both in our equation. We do that in \autoref{eq:basis sphere correction}. We view the energy reacing the atmosphere as a
circular area of energy, with the equation for the are of a circle being \autoref{eq:basis circle} \cite{areaCircle}. The area of a sphere is in \autoref{eq:basis sphere} \cite{areaSphere}. In
both equations, $r$ is the radius of the circle/sphere. By using \autoref{eq:basis circle} and \autoref{eq:basis sphere} in \autoref{eq:sb rewritten2} we get \autoref{eq:basis sphere2} where
$r$ is replaced by $R$. It is common in physics literature to use capitals for large objects like planets. However we are not done yet since we can divide some stuff out. We end up with
\autoref{eq:basis sphere final} as the final equation we are going to use.
\begin{subequations}
\begin{equation}
\pi r^2
\label{eq:basis circle}
\end{equation}
\begin{equation}
4 \pi r^2
\label{eq:basis sphere}
\end{equation}
\begin{equation}
\Delta T = \frac{\delta t (\pi R^2S - 4\pi R^2\sigma T^4)}{4\pi CR^2}
\label{eq:basis sphere2}
\end{equation}
\begin{equation}
\Delta T = \frac{\delta t (S - 4\sigma T^4)}{4C}
\label{eq:basis sphere final}
\end{equation}
\label{eq:basis sphere correction}
\end{subequations}
\subsection{Insolation}
With the current equation we calculate the global average surface temperature of the planet itself. However, this planet does not have an atmosphere just yet. Basically we modelled the
temperature of a rock floating in space, let's change that with \autoref{eq:atmos}. Here we assume that the area of the atmosphere is equal to the area of the planet surface. Obviously
this assumption is false, as the atmosphere is a sphere that is larger in radius than the planet, however the difference is not significant enough to account for it. We also define the
atmosphere as a single layer. This is due to the accessibility of the model, we want to make it accessible, not university simulation grade. One thing to take into account for the
atmosphere is that it only partially absorbs energy. The sun (or a similar star) is relatively hot and sends out energy waves (radiation) with relatively low wavelengths. The planet is
relatively cold and sends out energy at long wavelengths. As a side note, all objects radiate energy. You can verify this by leaving something in the sun on a hot day for a while and
almost touch it later. You can feel the heat radiating from the object. The planet is no exception and radiates heat as well, though at a different wavelength than the sun. The
atmosphere absorbs longer wavelengths better than short wavelengths. For simplicity's sake we say that all of the sun's energy does not get absorbed by the atmosphere. The planet's
radiation will be absorbed partially by the atmosphere. Some of the energy that the atmosphere absorbs is radiated into space and some of that energy is radiated back onto the planet's
surface. We need to adjust \autoref{eq:basis sphere final} to account for the energy being radiated from the atmosphere back at the planet surface.
So let us denote the temperature of the planet surface as $T_p$ and the temperature of the atmosphere as $T_a$. Let us also write the specific heat capacity of the planet surface as $C_p$
instead of $C$. We add the term in \autoref{eq:atmos on surface improved} to \autoref{eq:basis sphere final} in \autoref{eq:surface change}. In \autoref{eq:atmos on surface}, $\epsilon$ is the
absorbtivity of the atmosphere, the fraction of energy that the atmosphere absorbs. We divided \autoref{eq:atmos on surface} by $\pi R^2$ as we did that with \autoref{eq:basis sphere final} as
well, so we needed to make it match that division.
\begin{subequations}
\begin{equation}
4\pi R^2 \epsilon \sigma T_a^4
\label{eq:atmos on surface}
\end{equation}
\begin{equation}
4\epsilon \sigma T_a^4
\label{eq:atmos on surface improved}
\end{equation}
\begin{equation}
\Delta T_p = \frac{\delta t (S + 4\epsilon \sigma T_a^4 - 4\sigma T_p^4)}{4C_p}
\label{eq:surface change}
\end{equation}
\label{eq:atmos}
\end{subequations}
As you probably expected, the atmosphere can change in temperature as well. This is modelled by \autoref{eq:atmos change}, which is very similar to \autoref{eq:surface change}. There are
some key differences though. Instead of subtracting the radiated heat of the atmosphere once we do it twice. This is because the atmosphere radiates heat into space and towards the
surface of the planet, which are two outgoing streams of energy instead of one for the planet (as the planet obviously cannot radiate energy anywhere else than into the atmosphere).
$C_a$ is the specific heat capacity of the atmosphere.
\begin{equation}
\Delta T_a = \frac{\delta t (\sigma T_p^4 - 2\epsilon\sigma T_a^4)}{C_a}
\label{eq:atmos change}
\end{equation}
\subsection{The Latitude Longitude Grid}
With the current model, we only calculate the global average temperature. To calculate the temperature change along the surface and atmosphere at different points, we are going to use a grid.
Fortunately the world has already defined such a grid for us, the latitude longitude grid \cite{latlong}. The latitude is the coordinate running from the south pole to the north pole, with -90
being the south pole and 90 being the north pole. The longitude runs parallel to the equator and runs from 0 to 360 which is the amount of degrees that an angle can take when calculating the
angle of a circle. So 0 degrees longitude is the same place as 360 degrees longitude. To do this however we need to move on from mathematical formulae to code (or in this case pseudocode).
Pseudocode is a representation of real code. It is meant to be an abstraction of code such that it does not matter how you present it, but every coder should be able to read it and implement
it in their language of preference. This is usually easier to read than normal code, but more difficult to read than mathematical formulae. If you are unfamiliar with code or coding, look up
a tutorial online as there are numerous great ones.
The pseudocode in \autoref{alg:stream1v1} defines the main loop of the model. All values are initialised beforehand, based on either estimations, trial and error or because they are what they
are (like the Stefan-Boltzmann constant $\sigma$). The total time $t$ starts at 0 and increases by $\delta t$ after every update of the temperature. This is to account for the total time that
the model has simulated (and it is also used later). What you may notice is the $T_p[lan, lon]$ notation. This is to indicate that $T_p$ saves a value for each $lan$ and $lon$ combination.
It is initialised as all zeroes for each index pair, and the values is changed based on the calculations. You can view $T_p$ like the whole latitude longitude grid, where $T_p[lat, lon]$ is an
individual cell of that grid indexed by a specific latitude longitude combination.
\begin{algorithm}[hbt]
\SetAlgoLined
$\delta t \leftarrow 60 \cdot 5$ \;
$\sigma \leftarrow 5.67 \cdot 10^{-8}$ \;
$\epsilon \leftarrow 0.75$ \;
$C_p \leftarrow 10^7$ \;
$C_a \leftarrow 10^6$ \;
$S \leftarrow 1370$ \;
$R \leftarrow 6.4 \cdot 10^6$ \;
$t \leftarrow 0$ \;
\While{\texttt{TRUE}}{
\For{$lat \in [-90, 90]$}{
\For{$lon \in [0, 360]$}{
$T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t (S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p}$ \;
$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:stream1v1}
\end{algorithm}
\subsection{Day/Night Cycle}
As you can see, the amount of energy that reaches the atmopsphere is constant. However this varies based on the position of the sun relative to the planet. To fix this, we have to assign a function
to $S$ that gives the correct amount of energy that lands on that part of the planet surface. This is done in \autoref{alg:solar}. In this algorithm the term insolation is mentioned, which is $S$
used in the previous formulae if you recall. We use the $\cos$ function here to map the strength of the sun to a number between $0$ and $1$. The strength is dependent on the latitude, but since
that is in degrees and we need it in radians we transform it to radians by multiplying it by $\frac{\pi}{180}$. This function assumes the sun is at the equinox (center of the sun is directly
above the equator) \cite{equinox} at at all times. The second $\cos$ is needed to simulate the longitude that the sun has moved over the longitude of the equator. For that we need the difference
between the longitude of the point we want to calculate the energy for, and the longitude of the sun. The longitude of the sun is of course linked to the current time (as the sun is in a different\
position at 5:00 than at 15:00). So we need to map the current time in seconds to the interval $[0,$ seconds in a day$]$. Therefore we need the mod function. The mod function works like this:
$x$ mod $y$ means subtract all multiples of $y$ from $x$ such that $0 \leq x < y$. So to map the current time to a time within one day, we do $t$ mod $d$ where $t$ is the current time and $d$ is
the amount of seconds in a day. When we did the calculation specified in \autoref{alg:solar} we return the final value (which means that the function call is "replaced" \footnote{Replaced is not
necessarily the right word, it is more like a mathematical function $f(x)$ where $y = f(x)$. You give it an $x$ and the value that correpsonds to that $x$ is saved in $y$. So you can view the
function call in pseudocode as a value that is calculated by a different function which is then used like a regular number.} by the value that the function calculates). If the final value is less
than 0, we need to return 0 as the sun cannot suck energy out of the planet (that it does not radiate itself, which would happen if a negative value is returned).
In the second stream, it was revealed that $t$ mod $d$ in \autoref{alg:solar} should be $-t$ mod $d$ such that the sun moves in the right direction. In the first stream the sun would move to the
right (west to east), however the sun moves to the left (east to west) and so the time must be flipped in order for the model to be correct.
\begin{algorithm}[hbt]
\SetAlgoLined
\SetKwInput{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{insolation $ins$, latitude $lat$, longitude $lon$, time $t$, time in a day $d$}
\Output{Amount of energy $S$ that hits the planet surface at the given latitude-time combination.}
$longitude \leftarrow 360 \cdot \frac{(-t \text{ mod } d)}{d}$ \;
$S \leftarrow ins \cdot \cos(lat \frac{\pi}{180}) \cos((lon - longitude) \cdot \frac{\pi}{180})$ \;
\eIf{$S < 0$}{
\Return{$0$}
}{
\Return{$S$}
}
\caption{Calculating the energy from the sun (or similar star) that reaches a part of the planet surface at a given latitude and time}
\label{alg:solar}
\end{algorithm}
By implementing \autoref{alg:solar}, \autoref{alg:stream1v1} must be changed as well, as $S$ is no longer constant for the whole planet surface. So let us do that in \autoref{alg:stream1v2}. Note
that $S$ is defined as the call to \autoref{alg:solar} (as is showcased by the text \texttt{solar}). In case you are unfamiliar with calls, defining a function is defining how it works and
calling a function is actually using it.
\begin{algorithm}[hbt]
\SetAlgoLined
$\delta t \leftarrow 60 \cdot 5$ \;
$\sigma \leftarrow 5.67 \cdot 10^{-8}$ \;
$\epsilon \leftarrow 0.75$ \;
$C_p \leftarrow 10^7$ \;
$C_a \leftarrow 10^7$ \;
$I \leftarrow 1370$ \;
$R \leftarrow 6.4 \cdot 10^6$ \;
$t \leftarrow 0$ \;
$day \leftarrow 60 \cdot 60 \cdot 24$ \;
$S \leftarrow$ \texttt{solar($I$, $lat$, $lon$, $t$, $day$)} \;
$nlat$ is the amount of latitude points in the interval $[0, 90]$, how you divide them is your own choice. \;
$nlot$ is the amount of longitude points in the interval $[0, 360]$, how you divide them is your own choice. \;
\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 (S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p}$ \;
$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:stream1v2}
\end{algorithm}
\autoref{alg:stream1v2} calculates the values that are plotted (which is not discussed here as that is Python specific). Due to the \texttt{WHILE(TRUE)} loop, this calculation never finishes and
allows us to simulate days, weeks, months and even years of heat exchange all conveniently plotted in a graph. In Simon's implementation, the graphs update in realtime, meaning that whenever a
round of calculations has finished, they are immediately processed to be displayed in the graph.
However other forms of looking at the calculated data can be implemented, like writing a table to a txt file, saving the generated grpahs at a certain interval or spewing all the data into a csv
dataset. The possibilities are endless, and the whole goal of the model is for it to be modular. Meaning that if you want to do something with it (like have a multi-layered atmosphere instead of
a single layer atmosphere) you can just write some lines of code and run the model and it should still work. Therefore you can write your own extensions of the model to fit it to your needs and
requirements.

231
Documentation Source/Streams/Stream2.tex Normal file
View File

@ -0,0 +1,231 @@
\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}

266
Documentation Source/Streams/Stream3.tex Normal file
View File

@ -0,0 +1,266 @@
\section{Adding Mass to CLAuDE}
\subsection{The Momentum Equations}
The momentum equations are a set of equations that describe the flow of a fluid on the surface of a rotating body. For our model we will use the f-plane approximation. The equations corresponding
to the f-plane approximation are given in \autoref{eq:x momentum} and \autoref{eq:y momentum} \cite{momentumeqs}. Note that we are ignoring vertical moevement, as this does not have a significant
effect on the whole flow. All the symbols in \autoref{eq:x momentum} and \autoref{eq:y momentum} mean:
\begin{itemize}
\item $u$: The east to west velocity ($ms^{-1}$).
\item $t$: The time ($s$).
\item $f$: The coriolis parameter as in \autoref{eq:coriolis}.
\item $v$: The north to south velocity ($ms^{-1}$).
\item $\rho$: The density of the atmosphere ($kgm^{-3}$).
\item $p$: The atmospheric pressure ($Pa$).
\item $x$: The local longitude coordinate ($m$).
\item $y$: The local latitude coordinate ($m$).
\end{itemize}
If we then define a vector $\bar{u}$ as $(u, v, 0)$, we can rewrite both \autoref{eq:x momentum} as \autoref{eq:x momentum laplace}. Here $\nabla u$ is the gradient of
$u$ in both $x$ and $y$ directions. Then if we write out $\nabla u$ we get \autoref{eq:x momentum final}. Similarly, if we want to get $\delta v$ instead of $\delta u$ we rewrite
\autoref{eq:y momentum} to get \autoref{eq:y momentum laplace} and \autoref{eq:y momentum final}.
\begin{subequations}
\begin{equation}
\frac{Du}{Dt} - fv = -\frac{1}{\rho} \frac{\delta p}{\delta x}
\label{eq:x momentum}
\end{equation}
\begin{equation}
\frac{Dv}{Dt} - fu = -\frac{1}{\rho} \frac{\delta p}{\delta y}
\label{eq:y momentum}
\end{equation}
\begin{equation}
\frac{\delta u}{\delta t} + \bar{u} \cdot \nabla u - fv = -\frac{1}{\rho}\frac{\delta p}{\delta x}
\label{eq:x momentum laplace}
\end{equation}
\begin{equation}
\frac{\delta v}{\delta t} + \bar{u} \cdot \nabla v - fu = -\frac{1}{\rho}\frac{\delta p}{\delta y}
\label{eq:y momentum laplace}
\end{equation}
\begin{equation}
\frac{\delta u}{\delta t} + u\frac{\delta u}{\delta x} + v\frac{\delta u}{\delta y} - fv = -\frac{1}{\rho}\frac{\delta p}{\delta x}
\label{eq:x momentum final}
\end{equation}
\begin{equation}
\frac{\delta v}{\delta t} + u\frac{\delta v}{\delta x} + v\frac{\delta v}{\delta y} - fu = -\frac{1}{\rho}\frac{\delta p}{\delta y}
\label{eq:y momentum final}
\end{equation}
\end{subequations}
Now that we have the momentum equations sorted out, we need to define a method to do the gradient calculations for us. Therefore we define two functions \autoref{alg:gradient x} and
\autoref{alg:gradient y} that calculate the $x$ and $y$ gradients respectively.
\begin{algorithm}[hbt]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{Matrix (double array) $a$, first index $i$, second index $j$}
\Output{Gradient in the $x$ direction}
$grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon] - a[i, (j - 1) \text{ mod } nlon]}{\delta x[i]}$ \;
\Return{$grad$} \;
\caption{Calculating the gradient in the $x$ direction}
\label{alg:gradient x}
\end{algorithm}
\begin{algorithm}[hbt]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{Matrix (double array) $a$, first index $i$, second index $j$}
\Output{Gradient in the $y$ direction}
\eIf{$i == 0$ or $i == nlat - 1$}{
$grad \leftarrow 0$ \;
}{
$grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
}
\Return $grad$ \;
\caption{Calculating the gradient in the $y$ direction}
\label{alg:gradient y}
\end{algorithm}
With the gradient functions defined, we can move on to the main code for the momentum equations. The main loop is shown in \autoref{alg:stream3}. Do note that this loop replaces the one
in \autoref{alg:stream2v2} as these calculate the same thing, but the new algorithm does it better.
\begin{algorithm}
$S_{xu} \leftarrow \texttt{gradient\_x}(u, lan, lon)$ \;
$S_{yu} \leftarrow \texttt{gradient\_y}(u, lan, lon)$ \;
$S_{xv} \leftarrow \texttt{gradient\_x}(v, lan, lon)$ \;
$S_{yv} \leftarrow \texttt{gradient\_y}(v, lan, lon)$ \;
$S_{px} \leftarrow \texttt{gradient\_x}(p, lan, lon)$ \;
$S_{py} \leftarrow \texttt{gradient\_x}(p, lan, lon)$ \;
\While{\texttt{TRUE}}{
\For{$lat \in [1, nlat - 1]$}{
\For{$lon \in [0, nlon]$}{
$u[lan, lon] \leftarrow u[lan, lon] + \delta t \frac{-u[lan, lon]S_{xu} - v[lan, lon]S_{yu} + f[lan]v[lan, lon] - S_{px}}{\rho}$ \;
$v[lan, lon] \leftarrow v[lan, lon] + \delta t\frac{-u[lan, lon]S_{xv} - v[lan, lon]S_{yv} - f[lan]u[lan, lon] - S_{py}}{\rho}$ \;
}
}
}
\caption{Calculating the flow of the atmosphere (wind)}
\label{alg:stream3}
\end{algorithm}
\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). 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. $\Delta_x$ and $\Delta_y$ in \autoref{alg:laplacian} represents the calls to \autoref{alg:gradient x} and \autoref{alg:gradient y}
respectively. $\nabla^2$ in \autoref{alg:diffusion} represents the call to \autoref{alg:laplacian}.
\begin{algorithm}
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{A matrix (double array) a}
\Output{A matrix (double array) with results for the laplacian operator for each element}
\For{$lat \in [1, nlat - 1]$}{
\For{$lon \in [0, nlon]$}{
$output[lat, lon] \leftarrow \frac{\Delta_x(a, lat, (lon + 1) \text{ mod } nlon) - \Delta_x(a, lat, (lon - 1) \text{ mod } nlon)}{\delta x[lat]} + \frac{\Delta_y(a, lat + 1, lon) -
\Delta_y(a, lat - 1, lon)}{\delta y}$\;
}
}
\Return{$ouput$} \;
\caption{Calculate the laplacian operator over a matrix a}
\label{alg:laplacian}
\end{algorithm}
\begin{algorithm}
$\alpha_a \leftarrow 2 \cdot 10^{-5}$ \;
$\alpha_p \leftarrow 1.5 \cdot 10^{-6}$ \;
\While{\texttt{TRUE}}{
$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 loop for calculating the effects of diffusion}
\label{alg:diffusion}
\end{algorithm}
\subsection{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. 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.
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}
As we expect to use the divergence operator more often throughout our model, let us define a seperate function for it in \autoref{alg:divergence}. $\Delta_x$ and $\Delta_y$ in
\autoref{alg:divergence} represents the calls to \autoref{alg:gradient x} and \autoref{alg:gradient y} respectively. We do the multiplication with the velocity vector here already, as we expect
that we might use it in combination with the divergence operator more frequently.
\begin{algorithm}
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{A matrix (double array) $a$}
\Output{A matrix (double array) containing the result of the divergence operator taken over that element}
$dim_1 \leftarrow \text{ Length of } a \text{ in the first dimension}$ \;
\For{$i \in [0, dim_1]$}{
$dim_2 \leftarrow \text{ Length of } a \text{ in the second dimension (i.e. the length of the array stored at index } i)$ \;
\For{$j \in [0, dim_2]$}{
$output[i, j] \leftarrow \Delta_x(au, i, j) + \Delta_y(av, i, j)$ \;
}
}
\Return{$output$} \;
\caption{Calculate the result of the divergence operator on a vector}
\label{alg:divergence}
\end{algorithm}
With the divergence functon defined, 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}
$\alpha_a \leftarrow 2 \cdot 10^{-5}$ \;
$\alpha_p \leftarrow 1.5 \cdot 10^{-6}$ \;
\While{\texttt{TRUE}}{
$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 loop for calculating the effects of advection}
\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 chnage 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}
$\alpha_a \leftarrow 2 \cdot 10^{-5}$ \;
$\alpha_p \leftarrow 1.5 \cdot 10^{-6}$ \;
\While{\texttt{TRUE}}{
$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 loop for calculating the effects of advection}
\label{alg:advectionv2}
\end{algorithm}
Now that we have a varying density, we need to account for that in the temperature equations. So let us do that. We need it in the denominator as the density has a direct effect on the
heat capacity of the atmosphere. The changes are reflected in \autoref{alg:temperature with density}.
\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)}{\rho[lat, lon]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)}{\rho[lat, lon]C_a}$ \;
$t \leftarrow t + \delta t$ \;
}
}
}
\caption{The main loop of the temperature calculations}
\label{alg:temperature with density}
\end{algorithm}
\subsection{Improving the Coriolis Parameter}
Another change introduced is in the coriolis parameter. Up until now it has been a constant, however we know that it varies along the latitude. So let's make it vary over the latitude. Recall
\autoref{eq:coriolis}, where $\Theta$ is the latitude. Coriolis ($f$) is currently defined in \autoref{alg:gradient}, so let's incorporate the changes which are shown in \autoref{alg:coriolis}.
\begin{algorithm}
\SetAlgoLined
$C \leftarrow 2\pi R$ \;
$\delta y \leftarrow \frac{C}{nlat}$ \;
$\Omega \leftarrow 7.2921 \cdot 10^{-5}$ \;
\For{$lat \in [-nlat, nlat]$}{
$\delta x[lat] \leftarrow \delta y \cos(lat \cdot \frac{\pi}{180})$ \;
$f[lat] \leftarrow 2\Omega \sin(lat \frac{\pi}{180})$ \;
}
\caption{Calculating the gradient $\delta x$}
\label{alg:coriolis}
\end{algorithm}

83
Documentation Source/Streams/Stream4.tex Normal file
View File

@ -0,0 +1,83 @@
\section{Removing Some Assumptions and Mistakes from CLAuDE}
The first half of this stream was spent looking through the code and fixing some mistakes. To spare you dear reader from making these same mistakes, they have already been incorporated into
the previous sections, hooray! This does not only save you some work, but it also spares you from staring at a model that does not function due to wrongly defined constants or using the wrong
values.
\subsection{Adding a Spin-Up Time}
Instead of having a static start (having the planet start from rest, so no rotations allowed) we will have the model start up for some time before we start simulating the climate extensively.
To accomodate for this, we have to make some changes in the code. First we need to add two booleans (variables that can only take two values, either \texttt{TRUE} or \texttt{FALSE}) that we use
to indicate to the model whether we want to simulate the wind and whether we want to simulate advection. This means that the main loop will have some changes made to it. After performing the
calculations in \autoref{alg:temperature with density} we would calculate the velocities and afterwards we would calculate the advection. Instead let us change it to what is shown in
\autoref{alg:stream4v1}.
\begin{algorithm}
\While{\texttt{TRUE}}{
\autoref{alg:temperature with density} \;
\If{$velocity$}{
\autoref{alg:stream3} \;
\If{$advection$}{
\autoref{alg:advectionv2} \;
}
}
}
\caption{Main loop that can simulate flow and advection conditionally}
\label{alg:stream4v1}
\end{algorithm}
Now to dynamically enable/disable the simulation of flow and advection we need to add the spin-up calculations to the main loop. So in \autoref{alg:stream4v1}, before
\autoref{alg:temperature with density} we add \autoref{alg:spinup}. What it does is it changes the timestep when spinnning up and disables flow simulation, and when a week has passed it reduces
the timestep and enables flow simulation. At this point in time, the advection is not dynamically enabled/disabled but it is done by the programmer. Currently it will break the model, so I
recommend leaving it on \texttt{FALSE} until it is fixed in \autoref{sec:advectionfix}.
\begin{algorithm}
\eIf{$t < 7day$}{
$\delta t \leftarrow 60 \cdot 47$ \;
$velocity \leftarrow \texttt{FALSE}$ \;
}{
$\delta t \leftarrow 60 \cdot 9$ \;
$velocity \leftarrow \texttt{TRUE}$ \;
}
\caption{The spin-up block dynamically enabling or disabling flow simulation}
\label{alg:spinup}
\end{algorithm}
\subsection{Varying the Albedo}
The albdeo (reflectiveness of the planet's surface) is of course not the same over the whole planet. To account for this, we instead vary the albedo slightly for each point in the latitude
longitude grid. The algorithm that does this is shown in \autoref{alg:albedo variance}. The uniform distribution basically says that each allowed value in the interval has an equal chance of
being picked \cite{uniformdist}.
\begin{algorithm}
$V_a \leftarrow 0.02$ \;
\For{$lat \in [-nlat, nlat]$}{
\For{$lon \in [0, nlon]$}{
$R \leftarrow \text{ Pick a random number (from the uniform distribution) in the interval } [-V_a, V_a]$ \;
$a[lat, lon] \leftarrow a[lat, lon] + V_a \cdot R$\;
}
}
\caption{Varying the albedo of the planet}
\label{alg:albedo variance}
\end{algorithm}
\subsection{Fixing the Advection} \label{sec:advectionfix}
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}
$\alpha_a \leftarrow 2 \cdot 10^{-5}$ \;
$\alpha_p \leftarrow 1.5 \cdot 10^{-6}$ \;
$boundary \leftarrow 7$ \;
\While{\texttt{TRUE}}{
$T_{add} \leftarrow T_a + \delta t \alpha_a \nabla^2(T_a) + \nabla(T_a)$ \;
$T_a \leftarrow T_a - 0.5T_{add}[boundary:-boundary, :] \text{ //Only subtract } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + boundary, nlat - boundary]$. \;
$\rho[boundary: -boundary, :] \leftarrow \rho - 0.5(\delta t \nabla \rho) \text{ //Only change the density for indices in the interval } [-nlat + boundary, nlat - boundary]$ \;
$T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
}
\caption{The main loop for calculating the effects of advection}
\label{alg:advectionfix}
\end{algorithm}
\subsection{Adding Friction}
In order to simulate friction, we multiply the speeds $u$ and $v$ by $0.99$. Of course there are equations for friction but that gets complicated very fast, so instead we just assume that we
have a constant friction factor. This multiplication is done directly after \autoref{alg:stream3} in \autoref{alg:stream4v1}.

258
Documentation Source/Streams/Stream5.tex Normal file
View File

@ -0,0 +1,258 @@
\section{Up up and away! Adding More Layers to the Atmosphere}
Up until now we have neglected any vertical movement. This hampers the model, as the rising of warm air that then flows to the poles, cools down and flows back to the tropics is not possible
as the warm air cannot rise. So let us change that, let's add some vertical motion and some more layers to the atmosphere.
Remember \autoref{eq:atmos change}? We need this equation for every layer in the atmosphere. This also means that we have to adjust the main loop of the code, which is described in
\autoref{alg:temperature with density}. The $T_a$ needs to change, we need to either add a dimension (to indicate which layer of the atmosphere we are talking about) or we need to add different
matrices for each atmosphere layer. Let us define some useful variables in \autoref{alg:more layers}. We opt for adding a dimension as that costs less memory than defining new arrays
\footnote{This has to do with pointers, creating a new object always costs a bit more space than adding a dimension as we need a pointer to the object and what type of object it is whereas with
adding a dimension we do not need this additional information as it has already been defined}. So $T_a$, and all other matrices that have to do with the atmosphere (so not $T_p$ for instance)
are no longer indexed by $lat, lon$ but are indexed by $lat, lon, layer$.
\begin{algorithm}
$nlevels \leftarrow 4$ \;
$heights \leftarrow \text{Array with } nlevels \text{ layers, each with a uniform thickness of } \frac{100 \cdot 10^3}{nlevels} m$ \;
\caption{Definition of variables that are used throughout the code}
\label{alg:more layers}
\end{algorithm}
We also need to change all the gradient functions (\autoref{alg:gradient x} and \autoref{alg:gradient y}) to incorporate the atmospheric layers. Additionally we need a new gradient function that
calculates the gradient in the $z$ direction (vertical). Let us first change the existing gradient functions to take the atmospheric layer in effect. The changes can be found in
\autoref{alg:gradient x layer} and \autoref{alg:gradient y layer}. Let us improve the gradient in the $y$ direction as well. Since we are using the central difference method (calculating the
gradient by taking the difference of the next grid cell and the previous grid cell) there is no gradient at the poles. What we can do instead of returning $0$ for those cases is forward
differencing (calculating the gradient by taking the difference of the cell and the next/previous cell, multiplied by $2$ to keep it fair). A special change in both functions is checking whether
$k$ is equal to \texttt{NULL}. We do this as sometimes we want to use this function for matrices that does not have the layer dimension. Hence we define a default value for $k$ which is
\texttt{NULL}. \texttt{NULL} is a special value in computer science. It represents nothing. This can be useful sometimes if you declare a variable to be something but it is referring to something
that has been deleted or it is returned when some function fails. It usually indicates that something special is going on. So here we use it in the special case where we do not want to consider
the layer part in the gradient.
\begin{algorithm}[hbt]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$ with default value \texttt{NULL}}
\Output{Gradient in the $x$ direction}
\eIf{$k == \texttt{NULL}$}{
$grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon] - a[i, (j - 1) \text{ mod } nlon]}{\delta x[i]}$ \;
}{
$grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon, k] - a[i, (j - 1) \text{ mod } nlon, k]}{\delta x[i]}$ \;
}
\Return{$grad$} \;
\caption{Calculating the gradient in the $x$ direction}
\label{alg:gradient x layer}
\end{algorithm}
\begin{algorithm}[hbt]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$ with default value \texttt{NULL}}
\Output{Gradient in the $y$ direction}
\eIf{$k == \texttt{NULL}$}{
\uIf{$i == 0$}{
$grad \leftarrow 2 \frac{a[i + 1, j] - a[i, j]}{\delta y}$ \;
}\uElseIf{$i == nlat - 1$}{
$grad \leftarrow 2 \frac{a[i, j] - a[i - 1, j]}{\delta y}$ \;
}\uElse{
$grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
}
}{
\uIf{$i == 0$}{
$grad \leftarrow 2 \frac{a[i + 1, j, k] - a[i, j, k]}{\delta y}$ \;
}\uElseIf{$i == nlat - 1$}{
$grad \leftarrow 2 \frac{a[i, j, k] - a[i - 1, j, k]}{\delta y}$ \;
}\uElse{
$grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
}
}
\Return $grad$ \;
\caption{Calculating the gradient in the $y$ direction}
\label{alg:gradient y layer}
\end{algorithm}
With those changes done, let us define the gradient in the $z$ direction. The function can be found in \autoref{alg:gradient z layer}.
\begin{algorithm}[hbt]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$}
\Output{Gradient in the $z$ direction}
\uIf{$i == 0$}{
$grad \leftarrow 2\frac{a[i, j, k + 1] - a[i, j, k]}{\delta z[k]}$ \;
}\uElseIf{$i == nlat - 1$}{
$grad \leftarrow 2\frac{a[i, j, k] - a[i, j, k - 1]}{\delta z[k]}$ \;
}\uElse{
$grad \leftarrow \frac{a[i, j, k + 1] - a[i, j, k - 1]}{\delta z[k]}$ \;
}
\Return $grad$ \;
\caption{Calculating the gradient in the $z$ direction}
\label{alg:gradient z layer}
\end{algorithm}
As you can see, we have used $\delta z$ however, we have not defined it yet. Let us do that in \autoref{alg:gradient z}.
\begin{algorithm}[hbt]
\For{$k \in [0, nlevels - 1]$}{
$\delta z[k] \leftarrow heights[k + 1] - heights[k]$ \;
}
$\delta z[nlevels - 1] \leftarrow \delta z [nlevels - 2]$ \;
\caption{Defining $\delta z$ for later use throughout the code}
\label{alg:gradient z}
\end{algorithm}
Let's incorporate the changes for the Laplacian operator (\autoref{alg:laplacian}) as well. The new code can be found in \autoref{alg:laplacian layer}.
\begin{algorithm}[hbt]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{A matrix (double array) a}
\Output{A matrix (double array) with results for the laplacian operator for each element}
\eIf{$a.dimensions == 2$}{
\For{$lat \in [1, nlat - 1]$}{
\For{$lon \in [0, nlon]$}{
$output[lat, lon] \leftarrow \frac{\Delta_x(a, lat, (lon + 1) \text{ mod } nlon) - \Delta_x(a, lat, (lon - 1) \text{ mod } nlon)}{\delta x[lat]} + \frac{\Delta_y(a, lat + 1, lon) -
\Delta_y(a, lat - 1, lon)}{\delta y}$\;
}
}
}{
\For{$lat \in [1, nlat - 1]$}{
\For{$lon \in [0, nlon]$}{
\For{$k \in [0, nlevels - 1]$}{
$output[lat, lon, k] \leftarrow \frac{\Delta_x(a, lat, (lon + 1) \text{ mod } nlon, k) - \Delta_x(a, lat, (lon - 1) \text{ mod } nlon, k)}{\delta x[lat]} + \frac{\Delta_y(a,
lat + 1, lon, k) - \Delta_y(a, lat - 1, lon, k)}{\delta y} + \frac{\Delta_z(a, lat, lon, k + 1) - \Delta_z(a, lat, lon, k + 1)}{2\delta z[k]}$\;
}
}
}
}
\Return{$ouput$} \;
\caption{Calculate the laplacian operator over a matrix a}
\label{alg:laplacian layer}
\end{algorithm}
Of course we also need to incorporate the new layers in the divergence operator (\autoref{alg:divergence}). The new changes can be found in \autoref{alg:divergence layer}. Here we use $w$, the
vertical wind velocity. We define $w$ in the same way as $u$ and $v$, it is all zeroes (in the beginning) and has the same dimensions as $u$ and $v$.
\begin{algorithm}[!hbt]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{A matrix (double array) $a$}
\Output{A matrix (double array) containing the result of the divergence operator taken over that element}
$dim_1 \leftarrow \text{ Length of } a \text{ in the first dimension}$ \;
\For{$i \in [0, dim_1]$}{
$dim_2 \leftarrow \text{ Length of } a \text{ in the second dimension (i.e. the length of the array stored at index } i)$ \;
\For{$j \in [0, dim_2]$}{
$dim_3 \leftarrow \text{ Length of } a \text{ in the third dimension}$ \;
\For{$k \in [0, dim_3]$}{
$output[i, j] \leftarrow \Delta_x(au, i, j, k) + \Delta_y(av, i, j, k) + \Delta_z(aw, i, j, k)$ \;
}
}
}
\Return{$output$} \;
\caption{Calculate the result of the divergence operator on a vector}
\label{alg:divergence layer}
\end{algorithm}
With all those changes in the functions done, let us incorporate the changes into the model itself. We now need to account for the temperature change throughout the layers. Let us look at the
atmospheric temperature equation again (\autoref{eq:atmos change}). We need to account for one more thing, the absorbtion of energy from another layer. The new equation is shown in
\autoref{eq:atmos change layer}. Here $k$ is the layer of the atmosphere, $k = -1$ means that you use $T_p$ and $k = nlevels$ means that $T_{a_{nlevels}} = 0$ as that is space. Also, let us
rewrite the equation a bit such that the variables that are repeated are only written once and stuff that is divided out is removed, which is done in \autoref{eq:atmos change layer improved}.
Let us also clean up the equation for the change in the surface temperature (\autoref{eq:surface change}) in \autoref{eq:surface change improved}.
\begin{subequations}
\begin{equation}
\Delta T_{a_k} = \frac{\delta t (\sigma \epsilon_{k - 1}T_{a_{k - 1}}^4 + \sigma \epsilon_{k + 1}T_{a_{k + 1}}^4 - 2\epsilon_k\sigma T_{a_k}^4)}{C_a}
\label{eq:atmos change layer}
\end{equation}
\begin{equation}
\Delta T_{a_k} = \frac{\delta t \sigma (\epsilon_{k - 1}T_{a_{k - 1}}^4 + \epsilon_{k + 1}T_{a_{k + 1}}^4 - 2\epsilon_kT_{a_k}^4)}{C_a}
\label{eq:atmos change layer improved}
\end{equation}
\begin{equation}
\Delta T_p = \frac{\delta t (S + \sigma(4\epsilon_pT_a^4 - 4T_p^4))}{4C_p}
\label{eq:surface change improved}
\end{equation}
\end{subequations}
With the changes made to the equation, we need to make those changes in the code as well. We need to add the new dimension to all matrices except $T_p$ and $a$ as they are unaffected (with
regards to the storage of the values) by the addition of multiple atmospheric layers. Every other matrix is affected. The new code can be found in \autoref{alg:temperature layer}.
\begin{algorithm}[hbt]
\SetAlgoLined
\While{\texttt{TRUE}}{
\For{$lat \in [-nlat, nlat]$}{
\For{$lon \in [0, nlot]$}{
\For{$layer \in [0, nlevels]$}{
$T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t ((1 - a[lat, lon])S + \sigma(4\epsilon[0](T_a[lat, lon, 0])^4 - 4(T_p[lat, lon])^4))}
{4C_p[lat, lon]}$ \;
\uIf{$layer == 0$}{
$T_a[lat, lon, layer] \leftarrow T_a[lat, lon, layer] + \frac{\delta t \sigma((T_p[lat, lon])^4 - 2\epsilon[layer](T_a[lat, lon, layer])^4)}
{\rho[lat, lon, layer]C_a\delta z[layer]}$ \;
}\uElseIf{$layer == nlevels - 1$}{
$T_a[lat, lon, layer] \leftarrow T_a[lat, lon, layer] + \frac{\delta t \sigma(\epsilon[layer - 1](T_a[lat, lon, layer - 1])^4 - 2\epsilon[layer](T_a[lat, lon, layer])^4)}
{\rho[lat, lon, layer]C_a\delta z[layer]}$ \;
}\uElse{
$T_a[lat, lon, layer] \leftarrow T_a[lat, lon, layer] + \frac{\delta t \sigma(\epsilon[layer - 1](T_a[lat, lon, layer - 1])^4 + \epsilon[layer + 1]T_a[lat, lon, layer + 1]
- 2\epsilon[layer](T_a[lat, lon, layer])^4)}{\rho[lat, lon, layer]C_a\delta z[layer]}$ \;
}
$t \leftarrow t + \delta t$ \;
}
}
}
}
\caption{The main loop of the temperature calculations}
\label{alg:temperature layer}
\end{algorithm}
We also need to initialise the $\epsilon$ value for each layer. We do that in \autoref{alg:epsilon}.
\begin{algorithm}
$\epsilon[0] \leftarrow 0.75$ \;
\For{$i \in [1, nlevels]$}{
$\epsilon[i] \leftarrow 0.5\epsilon[i - 1]$
}
\caption{Intialisation of the insulation of each layer (also known as $\epsilon$)}
\label{alg:epsilon}
\end{algorithm}
Now we need to add vertical winds, or in other words add the $w$ component of the velocity vectors. We do that by editing \autoref{alg:stream3}. We change it to \autoref{alg:velocity}. Here we
use gravity ($g$) instead of the coriolis force ($f$) and calculate the pressure gradient in the $z$ direction.
\begin{algorithm}
$S_{xu} \leftarrow \texttt{gradient\_x}(u, lan, lon)$ \;
$S_{yu} \leftarrow \texttt{gradient\_y}(u, lan, lon)$ \;
$S_{xv} \leftarrow \texttt{gradient\_x}(v, lan, lon)$ \;
$S_{yv} \leftarrow \texttt{gradient\_y}(v, lan, lon)$ \;
$S_{px} \leftarrow \texttt{gradient\_x}(p, lan, lon)$ \;
$S_{py} \leftarrow \texttt{gradient\_y}(p, lan, lon)$ \;
$S_{pz} \leftarrow \texttt{gradient\_z}(p, lan, lon)$ \;
\While{\texttt{TRUE}}{
\For{$lat \in [1, nlat - 1]$}{
\For{$lon \in [0, nlon]$}{
\For{$layer \in [0, nlevels]$}{
$u[lan, lon, layer] \leftarrow u[lan, lon, layer] + \delta t \frac{-u[lan, lon, layer]S_{xu} - v[lan, lon, layer]S_{yu} + f[lan]v[lan, lon, layer] - S_{px}}{\rho}$ \;
$v[lan, lon, layer] \leftarrow v[lan, lon, layer] + \delta t \frac{-u[lan, lon, layer]S_{xv} - v[lan, lon, layer]S_{yv} - f[lan]u[lan, lon, layer] - S_{py}}{\rho}$ \;
$w[lan, lon, layer] \leftarrow w[lan, lon, layer] + \delta t (\frac{S_{pz}}{\rho[lan, lon, layer]} + g)$ \;
}
}
}
}
\caption{Calculating the flow of the atmosphere (wind)}
\label{alg:velocity}
\end{algorithm}
Lastly, 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
that the 3 dimensional matrix.
\begin{algorithm}
$\alpha_a \leftarrow 2 \cdot 10^{-5}$ \;
$\alpha_p \leftarrow 1.5 \cdot 10^{-6}$ \;
$boundary \leftarrow 7$ \;
\While{\texttt{TRUE}}{
$T_{add} \leftarrow T_a + \delta t \alpha_a \nabla^2(T_a) + \nabla(T_a)$ \;
$T_a \leftarrow T_a - 0.5T_{add}[boundary:-boundary, :, :] \text{ //Only subtract } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + boundary, nlat - boundary]$. \;
$\rho[boundary: -boundary, :, :] \leftarrow \rho - 0.5(\delta t \nabla \rho) \text{ //Only change the density for indices in the interval } [-nlat + boundary, nlat - boundary]$ \;
$T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
}
\caption{The main loop for calculating the effects of advection}
\label{alg:advection layer}
\end{algorithm}

24
Documentation Source/Streams/Stream6.tex Normal file
View File

@ -0,0 +1,24 @@
\section{Making a Dummy THICC Atmospheric Model*}
* The naming of this section is decided by the stream name, I did not come up with this \cite{twitch}. During this stream, a lot of plotting improvements have been made, which is not the scope of
this manual and hence has been left out. The plan was to add vertical momentum and advection, though things did not go according to plan\dots
\subsection{Discovering That Things Are Broken}
While trying to add vertical momentum, it appears that some parts of the model are broken in their current state. The horizontal advection is one of the things that is broken. If you recall,
we needed to use the Laplacian operator in the advection equations (as shown in \autoref{eq:diffusion}, diffusion is considered a part of advection since diffusion transports energy and matter
which is what advection does as well). The Laplacian operator (as shown in \autoref{alg:laplacian layer}) does not work. It is currently unknown why that is. If a fix has been found, the original
algorithm will be fixed accordingly (unless a lot of stuff changes and we need additional functions to be defined, then I will put a reference there to the section where it will be fixed). For
now we disabled the Laplacian operator as it has a small effect on the total advection (and because it is broken).
Another thing that we found out was broken is the vertical momentum. We tried to add it, ran into problems and just set it to 0 to fix the other problems that occured. One of those problems was
a wrong initialisation of the density. We basically told the model that the density is the same on every layer of the atmosphere, which is obviously not true. Hence we need to adjust that. The
new initialisation is described in \autoref{alg:density}. Note that the $\rho[:,: i]$ notation means that for every index in the first and second dimension, only change the value for the index $i$
in the third dimension.
\begin{algorithm}
$\rho[:, :, 0] \leftarrow 1.3$ \;
\For{$i \in [1, nlevels-1]$}{
$\rho[:, :, i] \leftarrow 0.1\rho[:, :, i - 1]$
}
\caption{Initialisation of the air density $\rho$}
\label{alg:density}
\end{algorithm}

47
Documentation Source/Streams/Stream7.tex Normal file
View File

@ -0,0 +1,47 @@
\section{Using Python to Model the Earth's Atmosphere}
This stream Simon was not feeling that well and it felt like his brain was not working, so be wary of errors! You have been warned. also the resolutin (size of an individual cell on the latitude
longitude grid) has been decreased to 5 degrees per cell instead of 3 degrees.
\subsection{Interpolating the Air Density}
In order to interpolate (see \autoref{sec:interpolation}) the air density, we need data. However currently we are just guessing the air density at higher levels, instead of taking real values.
So let us change that. For that we are going to use the U.S. Standard Atmosphere, an industry standard measure of the atmosphere on Earth \cite{usatmosp}. This data was provided in a text
(\texttt{TXT}) file which of course needs to be read in order for the data to be used in the model. Here we only care for the density and the temperature at a specific height. So the text file
only contains those two columns of the data (and the height in km of course as that is the index of the row, the property that uniquely identifies a row).
With that in mind, let's get coding and importing the data. We do this in \autoref{alg:usatmosp}. As one can see we do not specify how to open the file or how to split the read line, as this
is language specific and not interesting to describe in detail. I refer you to the internet to search for how to open a text file in the language you are working in. Keep in mind in which
magnitude you are working and in which magnitude the data is. If you work with $km$ for height and the data is in $m$, you need to account for that somewhere by either transforming the imported
data or work in the other magnitude.
\begin{algorithm}
$data \leftarrow \text{open text file containing the us standard atmosphere data}$ \;
\ForEach{$line \in data$}{
Split $line$ into three components, $sh, st$ and $sd$, representing the height, temperature and density respectively \;
$standardHeight.add(sh)$ \;
$standardTemperature.add(st)$ \;
$standardDensity.add(sd)$ \;
}
$densityProfile \leftarrow \texttt{interpolate}(heights, standardHeight, standardDensity)$ \;
$temperatureProfile \leftarrow \texttt{interpolate}(heights, standardHeight, standardTemperature)$ \;
\For{$alt \in [0, nlevels]$}{
$\rho[:, :, alt] \leftarrow densityProfile[alt]$ \;
$T_a[:, :, alt] \leftarrow temperatureProfile[alt]$ \;
}
\caption{Loading in the U.S. Standard Atmosphere}
\label{alg:usatmosp}
\end{algorithm}
Note that the function \texttt{interpolate} takes three arguments, the first one being the data points that we want to have values for, the second one is the data points that we know and the
third one is the values for the data points that we know. This function may or may not exist in your programming language of choice, which might mean that you have to write it yourself.
The formula that we use for interpolation can be found in \autoref{eq:interpolation}, though you still need to figure out what value you need for $\lambda$ (see \autoref{sec:interpolation}).
This is left as an exercise for the reader.
\subsection{Fixing Vertical Motion}
Another attempt was made at fixing the vertical motion. The changes are incorporated in \autoref{alg:advection layer}. Do keep in mind that the low air density in the upper layers messes a lot
with the vertical motion. In other words, it kinda works but not really. Another idea to help fix it, is to introduce a variable called $top$ which indicates the highest point that the
atmosphere may have. This value is initialised as $8 \cdot 10^3$ in meters (so 8 $km$). We then change the definition of $heights$ to: An array of uniform thickness of $\frac{top}{nlevels} m$.
We also added the $\delta z$ to \autoref{alg:temperature layer} as that was something that was still missing.
The current theory why the vertical velocity is not right is that the vertical thermodynamics may be wrong. This will be investigated further and we will report on this in future sections.

55
Documentation Source/Streams/TTNMETAF.tex Normal file
View File

@ -0,0 +1,55 @@
\appendix
\appendixpage
\section{Terms That Need More Explanation Then A Footnote}
\subsection{Potential} \label{sec:potential}
Potential is the energy change that occurs when the position of an object changes \cite{potential}. There are many potentials, like electric potential, gravitational potential and elastic
potential. Let me explain the concept with an example. Say you are walking on a set of stairs in the upwards direction. As your muscles move to bring you one step upwards, energy that is used
by your muscles is converted into gravitational potential. Now imagine you turn around and go downwards instead. Notice how that is easier? That is due to the gravitational potential being
converted back into energy so your muscles have to deliver less energy to get you down. The potential is usually tied to a force, like the gravitational force.
\subsection{Laplacian Operator} \label{sec:laplace}
The Laplacian operator ($\nabla^2$, sometimes also seen as $\Delta$) has two definitions, one for a vector field and one for a scalar field. The two concepts are not indpendent, a vector field
is composed of scalar fields \cite{vectorscalarfields}. Let us define a vector field first. A vector field is a function whose domain and range are a subset of the Eucledian $\mathbb{R}^3$ space.
A scalar field is then a function consisting out of several real variables (meaning that the variables can only take real numbers as valid values). So for instance the circle equation
$x^2 + y^2 = r^2$ is a scalar field as $x, y$ and $r$ are only allowed to take real numbers as their values.
With the vector and scalar fields defined, let us take a look at the Laplacian operator. For a scalar field $\phi$ the laplacian operator is defined as the divergence of the gradient of $\phi$
\cite{laplacian}. But what are the divergence and gradient? The gradient is defined in \autoref{eq:gradient} and the divergence is defined in \autoref{eq:divergence}. Here $\phi$ is a vector
with components $x, y, z$ and $\Phi$ is a vector field with components $x, y, z$. $\Phi_1, \Phi_2$ and $\Phi_3$ refer to the functions that result in the corresponding $x, y$ and $z$ values
\cite{vectorscalarfields}. Also, $i, j$ and $k$ are the basis vectors of $\mathbb{R^3}$, and the multiplication of each term with their basis vector results in $\Phi_1, \Phi_2$ and $\Phi_3$
respectively. If we then combine the two we get the Laplacian operator, as in \autoref{eq:laplacian scalar}.
\begin{subequations}
\begin{equation}
\text{grad } \phi = \nabla \phi = \frac{\delta \phi}{\delta x}i + \frac{\delta \phi}{\delta y}j + \frac{\delta \phi}{\delta z}k
\label{eq:gradient}
\end{equation}
\begin{equation}
\text{div} \Phi = \nabla \cdot \Phi = \frac{\delta \Phi_1}{\delta x} + \frac{\delta \Phi_2}{\delta y} + \frac{\delta \Phi_3}{\delta z}
\label{eq:divergence}
\end{equation}
\begin{equation}
\nabla^2 \phi = \nabla \cdot \nabla \phi = \frac{\delta^2 \phi}{\delta x^2} + \frac{\delta^2 \phi}{\delta y^2} + \frac{\delta^2 \phi}{\delta z^2}
\label{eq:laplacian scalar}
\end{equation}
\end{subequations}
For a vector field $\Phi$ the Laplacian operator is defined as in \autoref{eq:laplacian vector}. Which essential boils down to taking the Laplacian operator of each function and multiply it by
the basis vector.
\begin{equation}
\nabla^2 \Phi = (\nabla^2 \Phi_1)i + (\nabla^2 \Phi_2)j + (\nabla^2 \Phi_3)k
\label{eq:laplacian vector}
\end{equation}
\subsection{Interpolation} \label{sec:interpolation}
Interpolation is a form of estimation, where one has a set of data points and desires to know the values of other data points that are not in the original set of data points\cite{interpolation}.
Based on the original data points, it is estimated what the values of the new data points will be. There are various forms of interpolation like linear interpolation, polynomial interpolation
and spline interpolation. The CLAuDE model uses linear interpolation which is specified in \autoref{eq:interpolation}. Here $z$ is the point inbetween the known data points $x$ and $y$.
$\lambda$ is the factor that tells us how close $z$ is to $y$ in the interval $[0, 1]$. If $z$ is very close to $y$, $\lambda$ will have the value on the larger end of the interval, like 0.9.
Whereas if $z$ is close to $x$ then $\lambda$ will have a value on the lower end of the interval, like 0.1.
\begin{equation}
z = (1 - \lambda)x + \lambda y
\label{eq:interpolation}
\end{equation}

287
Documentation Source/references.bib Normal file
View File

@ -0,0 +1,287 @@
@misc {twitch,
howpublished="\url{https://www.twitch.tv/drsimonclark}",
journal={Twitch},
author={Clark, Simon}
}
@misc{playlist,
howpublished="\url{https://www.twitch.tv/collections/ewOcca6DExadGg}",
journal={Twitch},
author={Clark, Simon},
editor={Reed, CaptainEditor}
}
@misc{simon,
author={Clark, Simon}
}
@misc{SI,
title={SI Units},
howpublished="\url{https://www.nist.gov/pml/weights-and-measures/metric-si/si-units}",
journal={NIST},
publisher={National Institute of Standards and Technology, a U.S. Department of Commerce},
author={Isabel Chavez},
year={2019},
month={Nov}
}
%My Physics Book
@inbook{stefan-boltzmann,
place={Harlow},
title={Sears and Zemanskys University physics with modern physics},
publisher={Pearson Education},
author={Young, Hugh D. and Freedman, Roger A. and Ford, A. Lewis},
year={2016},
chapter={39},
pages={1328},
edition={14th global}
}
@inbook{thermo1,
place={Harlow},
title={Sears and Zemanskys University physics with modern physics},
publisher={Pearson Education},
author={Young, Hugh D. and Freedman, Roger A. and Ford, A. Lewis},
year={2016},
chapter={19},
pages={648},
edition={14th global}
}
@inbook{idealGas,
place={Harlow},
title={Sears and Zemanskys University physics with modern physics},
publisher={Pearson Education},
author={Young, Hugh D. and Freedman, Roger A. and Ford, A. Lewis},
year={2016},
chapter={18},
pages={610},
edition={14th global}
}
@inbook{densityAir,
place={Harlow},
title={Sears and Zemanskys University physics with modern physics},
publisher={Pearson Education},
author={Young, Hugh D. and Freedman, Roger A. and Ford, A. Lewis},
year={2016},
chapter={12},
pages={394},
edition={14th global}
}
@inbook{potential,
place={Harlow},
title={Sears and Zemanskys University physics with modern physics},
publisher={Pearson Education},
author={Young, Hugh D. and Freedman, Roger A. and Ford, A. Lewis},
year={2016},
chapter={7},
pages={227-247},
edition={14th global}
}
%Simon's physics book
@inbook{momentumeqs,
place={Cambridge},
edition={2nd},
DOI={10.1017/9781107588417.003},
title={Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation},
publisher={Cambridge University Press},
author={Vallis, Geoffrey K.},
year={2017},
chapter={2},
pages={69}
}
@inbook{diffusion,
place={Cambridge},
edition={2nd},
DOI={10.1017/9781107588417.003},
title={Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation},
publisher={Cambridge University Press},
author={Vallis, Geoffrey K.},
year={2017},
chapter={13},
pages={473}
}
@inbook{masscontinue,
place={Cambridge},
edition={2nd},
DOI={10.1017/9781107588417.003},
title={Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation},
publisher={Cambridge University Press},
author={Vallis, Geoffrey K.},
year={2017},
chapter={1},
pages={8}
}
@inbook{hydrostatic,
place={Cambridge},
edition={2nd},
DOI={10.1017/9781107588417.003},
title={Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation},
publisher={Cambridge University Press},
author={Vallis, Geoffrey K.},
year={2017},
chapter={1},
pages={12}
}
%Calculus
@inbook{areaCircle,
place={Don Mills, Ont.},
title={Calculus a complete course},
publisher={Pearson},
author={Adams, Robert Alexander and Essex, Christopher},
year={2018},
chapter={1},
pages={62},
edition={9th}
}
@inbook{areaSphere,
place={Don Mills, Ont.},
title={Calculus a complete course},
publisher={Pearson},
author={Adams, Robert Alexander and Essex, Christopher},
year={2018},
chapter={7},
pages={412},
edition={9th}
}
@inbook{vectorscalarfields,
place={Don Mills, Ont.},
title={Calculus a complete course},
publisher={Pearson},
author={Adams, Robert Alexander and Essex, Christopher},
year={2018},
chapter={15},
pages={867},
edition={9th}
}
@inbook{laplacian,
place={Don Mills, Ont.},
title={Calculus a complete course},
publisher={Pearson},
author={Adams, Robert Alexander and Essex, Christopher},
year={2018},
chapter={16},
pages={923},
edition={9th}
}
%Misc book sources
@inbook{uniformdist,
place={Hoboken, NJ},
edition={7th},
title={Applied statistics and probability for engineers},
publisher={Wiley},
author={Montgomery, Douglas C. and Runger, George C.},
year={2018},
chapter={3},
pages={49}
}
%General internet sources
@misc{latlong,
title={Geographic coordinate system},
howpublished="\url{https://en.wikipedia.org/wiki/Geographic_coordinate_system#Latitude_and_longitude}",
journal={Wikipedia},
publisher={Wikimedia Foundation},
author={Tobias Hoevekamp},
year={2020},
month={Jul}
}
@misc{equinox,
title={Equinox},
howpublished="\url{https://en.wikipedia.org/wiki/Equinox}",
journal={Wikipedia},
publisher={Wikimedia Foundation},
author={Palmen, Karl},
year={2020},
month={Jun}
}
@misc{mole,
title={SI Units - Amount of Substance},
howpublished="\url{https://www.nist.gov/pml/weights-and-measures/si-units-amount-substance}",
journal={NIST},
publisher={National Institute of Standards and Technology, a U.S. Department of Commerce},
author={Chavez, Isabel},
year={2019},
month={Nov}
}
@misc{specificGasConstantAir,
title={Gas constant},
howpublished="\url{https://en.wikipedia.org/wiki/Gas_constant#Specific_gas_constant}",
journal={Wikipedia},
publisher={Wikimedia Foundation},
author={192.2.69.128 (Anonymous Internet User)},
year={2020},
month={Jun}
}
@misc{coriolis,
title={Coriolis frequency},
howpublished="\url{https://en.wikipedia.org/wiki/Coriolis_frequency}",
journal={Wikipedia},
publisher={Wikimedia Foundation},
author={Kjkolb},
year={2020},
month={Jun}
}
@misc{circumference,
title={Circumference},
howpublished="\url{https://en.wikipedia.org/wiki/Circumference}",
journal={Wikipedia},
publisher={Wikimedia Foundation},
author={MBManie},
year={2020},
month={Jun}
}
@misc{albedo,
title={Albedo},
howpublished="\url{https://en.wikipedia.org/wiki/Albedo}",
journal={Wikipedia},
publisher={Wikimedia Foundation},
author={155.42.27.xxx (Anonymous Internet User)},
year={2020},
month={Jun}
}
@misc{fplane,
title={F-plane},
howpublished="\url{https://en.wikipedia.org/wiki/F-plane}",
journal={Wikipedia},
publisher={Wikimedia Foundation},
author={Johnson, Nathan},
year={2020},
month={Apr}
}
@book{usatmosp,
place={Washington, D.C.},
title={U.S. standard atmosphere},
publisher={National Oceanic and Atmospheric Administration},
author={National Aeronautics and Space Administration},
year={1976}
}
@misc{interpolation,
title={Interpolation},
howpublished = "\url{https://en.wikipedia.org/wiki/Interpolation}",
journal={Wikipedia},
publisher={Wikimedia Foundation},
author={Beldin, Dick},
year={2020},
month={May}
}