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

Merge pull request #1 from TechWizzart/rewrite 

Finished the rewrite so merge pls
This commit is contained in:
Sam Baggen 2020-09-08 20:20:21 +02:00 committed by GitHub
parent 37374321fe cf4c8b8311
commit c731b7256e
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
24 changed files with 1519 additions and 1542 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.

51
tex-docs/CLAuDE.tex
View File

@ -26,55 +26,54 @@
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}.
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}. There is a useful playlist on
Simon's Twitch which has all the streams without ad breaks or interruptions \cite{playlist}.
The manual itself is split up into distinct sections, each explaining one particular part of the model. Each section will be treating one topic, like radiation, advection or the control panel.
Although many concepts cannot be seen in isolation, as the wind has influence on how much temperature is distributed throughout the atmosphere, the calculations can be split up. The manual is
cumulative, starting with the basics and slowly building up to the current form of the algorithm. All changes to the algorithms can therefore be found here. An important distinction needs to be
made regarding the changes though. If the changes only change one part of the calculations, then it is considered an evolution, which will be added to the relevant section. However if the changes
are significant and not based on the previous code then the old alghorithms will be relocated to \autoref{sec:history}. Though the relevant theory will remain, as that is required to gain an
understanding of what the algorithm does. Do note that the radiation \autoref{sec:rad} is an exception for the first calculations as this forms the basis of the beginning of CLAuDE and the
fundamentals of the theory which I deem important enough to be left in place even if the calculations end up significantly different.
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
Within this manual we will not concern ourselves with plotting the data, instead we focus on the physics side of things and translating formulae 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
This manual is for the toy model, which is as of now still in development. One important thing to note is that 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
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.
Lastly, the manual is now up on the Planet Factory GitHub repository\cite{claudeGit}, together with all the source code. There is also a fork \cite{nomGit} that also contains the source code.
The manual is now up on the Planet Factory GitHub repository\cite{claudeGit}, together with all the source code. There is also a fork \cite{nomGit} that also contains the source code.
The fork will usually be more up to date than the version on the Planet Factory repository as Simon needs to merge pull requests into the repository. However I can update the fork freely so if a
particular stream is missing in the version on the Planet Factory repository, check the fork/Discord whether there is a newer version. If that is not the case, you just have to be a bit more
patient, or you can start writing a part of the manual yourself! Don't forget to ping me in the Discord to notify me of any additions (GitHub refuses to send me emails so I have no other way of
knowing).
\input{streams/Stream1.tex}
\input{topics/control_panel.tex}
\input{streams/Stream2.tex}
\input{topics/util_funcs.tex}
\input{streams/Stream3.tex}
\input{topics/radiation.tex}
\input{streams/Stream4.tex}
\input{topics/velocity.tex}
\input{streams/Stream5.tex}
\input{topics/advection.tex}
\input{streams/Stream6.tex}
\input{streams/Stream7.tex}
\input{streams/Stream8.tex}
\input{streams/Stream9.tex}
\input{streams/Stream10.tex}
\input{topics/master.tex}
\newpage
\input{streams/TTNMETAF.tex}
\input{appendices/TTNMETAF.tex}
\input{appendices/history.tex}
\input{appendices/vars.tex}
\newpage
\bibliography{references}

8
tex-docs/appendices/TTNMETAF.tex Normal file
View File

@ -0,0 +1,8 @@
\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.

138
tex-docs/appendices/history.tex Normal file
View File

@ -0,0 +1,138 @@
\section{History of the Algorithms} \label{sec:history}
Back when I was a young naive programmer, I made a thing. Now a few years down the line I made the thing again, but infinitely better. So I have no use for the old thing anymore. But fear not,
old algorithms (used by CLAuDE) will be collected here. This is just for historical purposes.
\subsection{Velocity}
\subsubsection{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
\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$ (note that this algorithm is obsolete)}
\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
\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).

54
tex-docs/appendices/vars.tex Normal file
View File

@ -0,0 +1,54 @@
\section{List of Variables}
Are you ever confused about what something is? Do you ever forget what a variable represents? Then I got the solution for you. The following overview will explain what each variable is and
represents. I will try to not use one variable for the same thing, though that is sometimes very difficult to do. I'll do my best. In the meantime, enjoy this exstensive list. Note that this
only applies to variables in code, every symbol in equations are explained at the equations themselves.
\begin{itemize}
\item $R$: The Gas Constant with value $8.3144621$ ($J(mol)^{-1}K$).
\item $day$: Length of one day in seconds ($s$).
\item $year$: Length of one year in seconds ($s$).
\item $\delta t$: How much time is between each calculation run in seconds ($s$).
\item $g:$ Magnitude of gravity on the planet in $ms^{-2}$.
\item $\alpha$: By how many degrees the planet is tilted with respect to the star's plane, also called axial tilt.
\item $top$: How high the top of the atmosphere is with respect to the planet surface in meters ($m$).
\item $ins$: Amount of energy from the star that reaches the planet per unit area ($Jm^{-2}$).
\item $\epsilon$: Absorbtivity of the atmosphere, fraction of how much of the total energy is absorbed (unitless).
\item $resolution$: The amount of degrees on the latitude longitude grid that each cell has, with this setting each cell is 3 degrees latitude high and 3 degrees longitude wide.
\item $nlevels \leftarrow 10$: The amount of layers in the atmosphere.
\item $\delta t_s$: The time between calculation rounds during the spin up period in seconds ($s$).
\item $t_s$: How long we let the planet spin up in seconds ($s$).
\item $adv$: Whether we want to enable advection or not.
\item $velocity$: Whether we want to calculate the air velocity or not.
\item $adv\_boun$: How many cells away from the poles where we want to stop calculating the effects of advection.
\item $nlon$: The amount of longitude gridpoints that we use, which depends on the resolution.
\item $nlat$: The amount of latitude gridpoints that we use, which depends on the resolution.
\item $T_p$: The temperature of the planet, a 2D array representing a latitude, longitude grid cell.
\item $T_a$: The temperature of the atmosphere, a 3D array representing a grid cell on the latitude, longitude, atmospheric layer grid.
\item $\sigma$: The Stefan-Boltzmann constant equal to $5.670373 \cdot 10^{-8} \ (Wm^{-2}K^{-4})$.
\item $C_a$: Specific heat capacity of the air, equal to $1.0035 Jg^{-1}K^{-1}$.
\item $C_p$: Specific heat capacity of the planet, equal to $1.0 \cdot 10^{6} Jg^{-1}K^{-1}$.
\item $a$: Albedo, the reflectiveness of a substance. Note that $a$ is used in general functions as an array that is supplied as input. If that is the case it can be read at the top of the
algorithm.
\item $\rho$: The density of the atmosphere, a 3D array representing a grid cell on the latitude, longitude, atmospheric layer grid.
\item $\delta x$: How far apart the gridpoints are in the $x$ direction in degrees longitude.
\item $\delta y$: How far apart the gridpoints are in the $y$ direction in degrees latitude.
\item $\delta z$: How far apart the gridpoints are in the $z$ direction in $m$.
\item $heights$: How high an atmospheric layer is in $m$.
\item $\tau$: The optical depth for an atmospheric layer.
\item $\tau_0$: The optical depth at the planet surface.
\item $f_l$: The optical depth parameter.
\item $pressureProfile$: The average pressure taken over all atmospheric layers in a latitude, longitude gridcell.
\item $densityProfile$:The average density taken over all atmospheric layers in a latitude, longitude gridcell.
\item $temperatureProfile$: The average temperature taken over all atmospheric layers in a latitude, longitude gridcell.
\item $U$: Upward flux of radiation, 1D array representing an atmospheric layer.
\item $D$: Downward flux of radiation, 1D array representing an atmospheric layer.
\item $u$: The east to west air velocity in $ms^{-1}$.
\item $v$: The north to south air velocity in $ms^{-1}$.
\item $w$: The bottom to top air velocity in $ms^{-1}$.
\item $f$: The coriolis parameter.
\item $\Omega$: The rotation rate of the planet in $rads^{-1}$.
\item $p$: The pressure of a latitude, longitude, atmospheric layer gridcell.
\item $p_0$: The pressure of a latitude, longitude, atmospheric layer gridcell from the previous calculation round.
\item $\alpha_a$: The thermal diffusivity constant for air.
\item $\alpha_p$: The thermal diffusivity constant for the planet surface.
\end{itemize}

11
tex-docs/references.bib
View File

@ -45,6 +45,17 @@ pages={1328},
edition={14th global}
}
@inbook{specificHeat,
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={17},
pages={581},
edition={14th global}
}
@inbook{thermo1,
place={Harlow},
title={Sears and Zemanskys University physics with modern physics},

29
tex-docs/streams/Stream.aux
View File

@ -1,29 +0,0 @@
\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
tex-docs/streams/Stream1.tex
View File

@ -1,228 +0,0 @@
\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.

81
tex-docs/streams/Stream10.tex
View File

@ -1,81 +0,0 @@
\section{Putting Our Homemade Climate Model Through Its Paces}
Big stream this stream as we got some working code! Always great when stuff works. This stream, we tackled the radiation problem, added axial tilt to the planet, fixed vertical motion (but not
advection), added stratospheric heating and some other code clean up stuff. This means that the big rework is getting closer and closer. How exciting!
%TESTSSSS WE GOT WORKING CODE! HELLO WORLD!
\subsection{Fixing Up the Code}
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 differenc 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$ funciton in \autoref{alg:divergence layer}. 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.
We also fixed up the radiation scheme, as shown in \autoref{alg:optical depth}. Basically we had the definition of $U[k + 1] = \text{something something} U[k + 1]$. This means that the definition
was relying on itself, which is obviously impossible and wrong. So we changed it to it's current form and it is fixed, hooray!
Vertical motion has also been fixed, as shown in \autoref{alg:velocity}. Due to some error in the representation of the vertical motion it did not work. So we changed from that representation to
another. Now the vertical velocity is proportional to the rate of change of the pressure which does work like it should.
\subsection{Tilting the Planet}
In order to model a planet that has seasons, like Earth, we need to tilt the planet. This has as effect that the sun is not always directly above the equator but is above a certain band around
the equator as the year moves on. This means that some hemispheres receive more/less sun based on what part of the year it is. Which corresponds to the various seasons we have on Earth. But in
order to do that, we have to change the \texttt{solar} function. The new version as shown in \autoref{alg:solar tilt} will replace \autoref{alg:solar}. Here $\alpha$ is the tilt in degrees.
\begin{algorithm}
\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, longitude and time combination.}
$sun\_lon \leftarrow -t \text{ mod } d$ \;
$sun\_lon \leftarrow sun\_lon \cdot \frac{360}{d}$ \;
$sun\_lat \leftarrow \alpha\cos(\frac{2t\pi}{year})$ \;
$S \leftarrow insolation\cos(\frac{\pi(lat - sun\_lat)}{180})$ \;
\uIf{$S < 0$}{
\Return{$0$} \;
} \uElse {
$lon\_diff \leftarrow lon - sun\_lon$ \;
$S \leftarrow S\cos(\frac{lon\_diff\pi}{180})$ \;
\uIf{$S < 0$}{
\uIf{$lat + sun\_lat > 90$ or $lat + sun\_lat < -90$}{
\Return{$insolation\cos(\frac{\pi(lat + sun\_lat)}{180})\cos(\frac{lon\_diff\pi}{180})$} \;
} \uElse {
\Return{$0$} \;
}
} \uElse {
\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 tilt}
\end{algorithm}
What the code in \autoref{alg:solar tilt} does boils down to calculating the latitude and longitude of the sun and checking whether the planet receives any energy. If not return $0$ immediately.
If so we check if the difference between the sun's longitude and the planet's longitude and calculate how much energy would hit the planet given that the sun is not straight above the equator.
We do this by multiplying the energy it would receive from the sun if it were above the equator $S$ by the cosine of the difference in longitudes, which represents the tilt. Then we check again
if the planet is receiving energy, if not we check if it happens around the poles. We do this because due to the tilt it can be the case that at certain points in the year the pole is in constant
sunlight, i.e. the sun does not go down. This creates a sort of overshoot which needs to be accounted for. If it does this then we add the latitudes of the sun and the planet together and use
that to calculate the energy that would hit that spot. If it is not the case that we are around the poles and we do not receive energy, then we return $0$. If it happens to be that we do receive
energy (so no negative values) then we return $S$.
\subsection{Adding In Some Ozone (Or Something Else That Approximates It)}
Adding in ozone in the stratosphere is hella complicated, so we leave that as an exercise to the reader as in true academic fashion. Just joking, if you want you can work on implementing ozone
however we opt not to because it is quite complicated. Instead we approximate it, which is decent enough for our purpose. We need to do it in \autoref{alg:optical depth} as we need to adjust the
$Q$. We add in a check to see if we are currently calculating the radiation in the stratosphere. If so we add some radiation extra to replicate the effect of ozone. As can be seen in
\autoref{alg:ozone}, where we only focus on the $Q$ part of \autoref{alg:optical depth}, we add in some extra radiation based on how high the current layer calculation is, which scales with the
height.
\begin{algorithm}
\For{$level \in [0, nlevels]$}{
$Q[level] \leftarrow - \frac{S_z(U - D, 0, 0, level)}{10^3 \cdot densityProfile[level]}$ \;
\uIf{$heights[level] > 20 \cdot 10^3$}{
$Q[level] \leftarrow Q[level] + \texttt{solar}(5, lat, lon, t) \frac{24 \cdot 60 \cdot 60(\frac{heights[level] - 20 \cdot 10^3}{10^3})^2}{30^2}$ \;
}
}
\caption{Replicating the effect of ozone}
\label{alg:ozone}
\end{algorithm}
It is at this point that we reached the state that CLAuDE is in a testable state. This means that we have the model working in such a way that we can do some simple experiments like altering how
long a day is, what would happen if the sun would send out more energy (which usually means that it is bigger) or what would happen if you tidally lock a planet (stop it rotating completely).

231
tex-docs/streams/Stream2.tex
View File

@ -1,231 +0,0 @@
\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
tex-docs/streams/Stream3.tex
View File

@ -1,266 +0,0 @@
\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 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}
$\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
tex-docs/streams/Stream4.tex
View File

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

275
tex-docs/streams/Stream5.tex
View File

@ -1,275 +0,0 @@
\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}. Here $a.dimensions$ is the attribute that tells us how
deeply nested the array $a$ is. If the result is $1$ we have just a normal array, if it is $2$ we have a double array (an array at each index of the array) which is also called a matrix and if it
is $3$ we have a triple array. We need this because we have a one-dimensional case, for when we do not use multiple layers and a three-dimensional case for when we do use multiple layers. This
distinction is needed to avoid errors being thrown when running the model with one or multiple layers.
\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{$a.dimensions == 1$}{
\uIf{$k == 0$}{
$grad \leftarrow \frac{a[k + 1] - a[k]}{\delta z[k]}$ \;
}\uElseIf{$k == nlevels - 1$}{
$grad \leftarrow \frac{a[k] - a[k - 1]}{\delta z[k]}$ \;
}\uElse{
$grad \leftarrow \frac{a[k + 1] - a[k - 1]}{2\delta z[k]}$ \;
}
} \uElse {
\uIf{$k == 0$}{
$grad \leftarrow \frac{a[i, j, k + 1] - a[i, j, k]}{\delta z[k]}$ \;
}\uElseIf{$k == nlevels - 1$}{
$grad \leftarrow \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]}{2\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 change in pressure. Therefore we need to store a copy of the pressure before we do any calculations. This needs to be a
copy due to aliasing \footnote{Aliasing is assigning a different name to a variable, while it remains the same variable. Take for instance that we declare a variable $x$ and set it to be $4$.
Then we say $y \leftarrow x$, which you might think is the same as saying they $y \leftarrow 4$ but behind the screen it is pointing to $x$. So if $x$ changes, then so does $y$.}
\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)$ \;
\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[lat, lon, layer] + \delta t \frac{-u[lat, lon, layer]S_{xu} - v[lat, lon, layer]S_{yu} + f[lat]v[lat, lon, layer] - S_{px}}{\rho}$ \;
$v[lan, lon, layer] \leftarrow v[lat, lon, layer] + \delta t \frac{-u[lat, lon, layer]S_{xv} - v[lat, lon, layer]S_{yv} - f[lat]u[lat, lon, layer] - S_{py}}{\rho}$ \;
$w[lan, lon, layer] \leftarrow w[lat, lon, layer] - \frac{p[lat, lon, layer] - p_o[lat, lon, layer]}{\delta t\rho[lat, lon, layer]g}$ \;
}
}
}
$p_o \leftarrow copy(p)$ \;
}
\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
of 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}

25
tex-docs/streams/Stream6.tex
View File

@ -1,25 +0,0 @@
\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}) did not work. This is because there was a misplaced bracket causing weird numerical errors.
This has been fixed in the code (but was never present in the manual, yay for me) and can be safely enabled, though for this stream we disabled the Laplacian operator as it has a small effect on
the total advection (and because it was at this time broken). %The Laplacian has been re-enabled in Stream 8 (\autoref{sec:stream8}).
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. As soon as the vertical momentum has been fixed it will be fixed in the correct spot. If a lot of the code changes then it will probably be in another section and I will
insert a reference to that section.
\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
tex-docs/streams/Stream7.tex
View File

@ -1,47 +0,0 @@
\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.

117
tex-docs/streams/Stream8.tex
View File

@ -1,117 +0,0 @@
\section{Getting Radiation Right in our Climate Model! 3D Motion Here We Come} \label{sec:stream8}
The time has come to finally fix 3D motion. For this to work, we need to use a radiation scheme, which Simon \sout{shamelessly stole} got inspired by the Isca project \cite{isca}. So he followed
the references and found a paper which he is going to use in our model \cite{greyRad}. Great, so let's get into it shall we.
\subsection{Grey Radiation Scheme}
A radiation scheme is a model for how energy is redistributed using light in a system. Such a model is a Grey radiation scheme if you split it into two parts, short and long wavelength radiation.
So you have two redistribution systems, one for short wavelength light and one for long wavelength light. Another assumption we make when using the Grey radiation scheme, is that the atmosphere
is transparent to short wavelength radiation, meaning it lets through light with short wavelengths. Additionally we use a two stream approximation, which means that we have a stream of radiation
going up, and another stream of radiation going down. Note that these two streams are both long wavelength radiation, because we said earlier we assume the atmosphere completely ignores short
wavelength radiation.
The two long wavelength radiation streams are described in \autoref{eq:upward radiation} and \autoref{eq:downward radiation} \cite{greyRad}. In those equations, the symbols are:
\begin{itemize}
\item $U$: Upward flux.
\item $D$: Downward flux.
\item $B$: The Stefan-Boltzmann equation (see \autoref{eq:stefan-boltzmann}).
\item $\tau$: Optical depth.
\end{itemize}
\begin{subequations}
\begin{equation}
\frac{dU}{d\tau} = U - B
\label{eq:upward radiation}
\end{equation}
\begin{equation}
\frac{dD}{d\tau} = B - D
\label{eq:downward radiation}
\end{equation}
\end{subequations}
With \autoref{eq:upward radiation} and \autoref{eq:downward radiation} written down, we can discuss how they work. These equations need a boundary condition to work, a starting point if you like.
For those equations the boundary conditions are that $U$ is at the surface equal to $B$ and that $D$ at the top of the atmosphere is equal to $0$. Meaning that in the beginning the top of the
atmosphere has no downward flux as there is no heat there, and that the bottom of the atmosphere has a lot of upward flux as most if not all of the heat is located there. Then after the spin up
time this should stabilise. We are interested in the change of the fluxes, so $dU$ and $dD$, to get those we need to multiply the right hand side by $d\tau$. Then we have the flow of radiation
that we need. However we cannot solely use these two equations to calculate the heat of a given layer. For that we need a few more components. These are described in \autoref{eq:heat layer}.
Here $Q_R$ is the amount of heat in a layer, $c_p$ is the specific heat capacity of dry air (our atmosphere), $\rho$ is the density of the air in that layer and $\delta z$ is the change in height.
$\delta U - D$ are the change in net radiation, meaning the amount of radiation that is left over after you transferred the upward and downward flux. See it as incoming and outgoing energy for a
given layer, the net change (either cooling down or heating up) is what remains after you have subtracted the incoming energy from the outgoing energy. While this explanation is not entirely true
(as flux is not entirely equivalent to energy), it explains the concept the best.
\begin{equation}
Q_R = \frac{1}{c_p\rho}\frac{\delta(U - D)}{\delta z}
\label{eq:heat layer}
\end{equation}
Now only one question remains: what is optical depth? Optical depth is the amount of work a photon has had to do to get to a certain point. This might sound really vague, but bear with me.
Optical depth describes how much stuff a certain photon has had to go through to get to a point. As you'd expect this is $0$ at the top of the atmosphere as space is a big vacuum so no stuff to
move through, so no work. Then the further the photon moves into the atmosphere, the more work the photon has had to do to get there. This is because it now needs to move through gases, like air,
water vapour and other gases. Hence the closer the photon gets to the surface of the planet, the larger the optical depth is because the photon has had to work more to get there. This phenomenon
is described in \autoref{eq:optical depth}. The symbols in the equation mean:
\begin{itemize}
\item $\tau_0$: Optical depth at the surface.
\item $p$: Atmospheric pressure ($Pa$).
\item $p_s$: Atmospheric pressure at the surface ($Pa$).
\item $f_l$: The linear optical depth parameter, with a value of 0.1.
\end{itemize}
\begin{equation}
\tau = \tau_0[f_l(\frac{p}{p_s}) + (1 - f_l)(\frac{p}{p_s})^4]
\label{eq:optical depth}
\end{equation}
As one can see, \autoref{eq:optical depth} has two parts, a linear part and a quatric part (to the power $4$). The quatric term approximates the structure of water vapour in the atmosphere, which
roughly scales with $\frac{1}{4}$ with respect to the height. The linear term is present to fix numerical behaviour because this is an approximation which will not be completely correct (that's
why it is an approximation) so we add this term to make it roughly right. The same thing holds for $f_l$ which can be manually tuned to fix weird numerical behaviour.
\subsection{Getting the equations to code}
With these equations in our mind, let's get coding. First we add the pressure profile, the pressure of all atmospheric layers at a lat lon point. We need this to accurately represent the optical
depth per atmospheric layer. Then we need to use the pressure profile with regards to \autoref{eq:optical depth}. The resulting code can be found in \autoref{alg:optical depth}. This algorithm
replaces the temperature calculations we have done in \autoref{alg:temperature layer}, as this is basically a better version of the calculations done in that algorithm. $f_l$ has a value of $0.1$
and is located near all the other constants in the code, henceforth we will refer to this section in the code as the control panel, since most if not all of the constants can be tweaked here.
$\tau_0$ is a function that gives the surface optical depth for a given latitude. The corresponding equation can be found in \autoref{eq:optical depth surface} \cite{simon}. Translating this
into code is left as an exercise to the reader. $U[0]$ is the boundary condition discussed before (being the same as \autoref{eq:stefan-boltzmann}), just as $D[nlevels]$ is the boundary condition.
$S_z$ represents the call to \autoref{alg:gradient z layer}. \texttt{solar} represents the call to \autoref{alg:solar}.
\begin{algorithm}
\For{$lat \in [-nlat, nlat]$}{
\For{$lon \in [0, nlon]$}{
$pressureProfile \leftarrow p[i,j,:]$ \;
$\tau = \tau_0(lat)f_l\frac{pressureProfile}{pressureProfile[0]} + (1 - f_l)(\frac{pressureProfile}{pressureProfile[0]})^4)$ \;
$U[0] \leftarrow \sigma T_p[lat, lon]^4$ \;
\For{$level \in [1, nlevels]$}{
$U[level] \leftarrow U[level - 1] - \frac{(\tau[level] - \tau[level - 1])(\sigma \cdot (mean(T_a[:, :, level]))^4)}{1 + (\tau[level - 1] - \tau[level])}$ \;
}
$D[nlevels - 1] \leftarrow 0$ \;
\For{$level \in [nlevels - 1, 0]$}{
$D[level] \leftarrow D[level + 1] - \frac{(\tau[level + 1] - \tau[level])(\sigma \cdot (mean(T_a[:, :, level]))^4)}{1 + (\tau[level] - \tau[level + 1])}$ \;
}
\For{$level \in [0, nlevels]$}{
$Q[level] \leftarrow - \frac{S_z(U - D, 0, 0, level)}{10^3 \cdot densityProfile[level]}$ \;
}
$T_a[lat, lon, :] \leftarrow T_a[lat, lon, :] + Q$ \;
$S \leftarrow \texttt{solar}(I, lat, lon, t)$ \;
$T_p[lat, lon] \leftarrow T_p[lat, lon] \frac{\delta t((1 - a[lat, lon]) S + S_z(D, 0, 0, 0) - \sigma T_p[lat, lon]^4)}{C_p[lat ,lon]}$ \;
}
}
\caption{Adding in radiation}
\label{alg:optical depth}
\end{algorithm}
\begin{equation}
\tau_0 = 3.75 + \cos(lat \frac{\pi}{90})\frac{4.5}{2}
\label{eq:optical depth surface}
\end{equation}
%Note that in this form, it did not work on stream yet. This may be due to a coding error or to a missing equation, constant or something similar. If it turns out to be a simple fix, then it will
%be fixed in this section. If a lot of other things change in order for the fix to work, then it will probably be a seperate section with a reference to that section here.

79
tex-docs/streams/Stream9.tex
View File

@ -1,79 +0,0 @@
\section{Starting to Deal With the Poles}
It is time to deal with the pole situation. The north and south poles that is, not the lovely people over in Poland. We run into problems because the latitude longitude grid cells become to small
near the poles. Therefore, the magnitudes no longer fit into one cell and overflow into other cells which makes everything kind of funky. So we need to fix that, and we do that by a planar
approximation.
\subsection{The Theory Behind the Planar Approximation}
As said earlier, the grid cells on the latitude longitude grid get closer together the closer you get to the poles which poses problems. To fix this, we will be using a planar approximation of
the poles. What this means is that we will map the 3D grid near the poles onto a 2D plane parallel to the poles, as if we put a giant flat plane in the exact center of the poles and draw lines
from the grid directly upwards to the plane. For a visual representation, please consult the stream with timestamp 1:38:25 \cite{polarPlane}, which includes some explanation. In the streamm we
use $r$ to indicate the radius of the planet (which we assume is a sphere), $\theta$ for the longitude and $\lambda$ for the latitude. So we have spherical coordinates, which we need to transform
into $x$ and $y$ coordinates on the plane. We also need the distance between the center point (the point where the plane touches the planet which is the center of the pole) and the projected
point on the plane from the grid (the location on the plane where a line from the gird upwards to the plane hits it). This distance is denoted by $a$ (Simon chose this one, not me). We then get
the following equations as shown in \autoref{eq:polar distance}, \autoref{eq:polar x} and \autoref{polar y}.
\begin{subequations}
\begin{equation}
a = r \cos(\theta)
\label{eq:polar distance}
\end{equation}
\begin{equation}
x = a \sin(\lambda)
\label{eq:polar x}
\end{equation}
\begin{equation}
y = a \cos(\lambda)
\label{eq:polar y}
\end{equation}
\end{subequations}
But what if we know $x$ and $y$ and want to know $\theta$ and $\lambda$? Pythagoras' Theorem then comes into play \cite{pythagoras}. We know that (due to Pythagoras) \autoref{eq:pythagoras} must
always be true. Then if we substitue $a$ by $\sqrt{x^2 + y^2}$ in \autoref{eq:polar distance} we get \autoref{eq:polar theta1}. Then we transform that equation such that we only have $\theta$ on
one side and the rest on the other side (since we want to know $\theta$) and we get \autoref{eq:polar theta3}.
\begin{equation}
x^2 + y^2 = a^2
\label{eq:pythagoras}
\end{equation}
\begin{subequations}
\begin{equation}
\sqrt{x^2 + y^2} = r\cos(\theta)
\label{eq:polar theta1}
\end{equation}
\begin{equation}
\frac{\sqrt{x^2 + y^2}}{r} = \cos(\theta)
\label{eq:polar theta2}
\end{equation}
\begin{equation}
\cos^{-1}(\frac{\sqrt{x^2 + y^2}}{r}) = \theta
\label{eq:polar theta3}
\end{equation}
\end{subequations}
For $\lambda$ we need another trigonometric function which is the tangent ($\tan$). The tangent is defined in \autoref{eq:tan}. If we then take a look at \autoref{eq:polar x} and
\autoref{eq:polar y}, we see that $\lambda$ is present in both equations. So we need to use both to get $\lambda$ \footnote{Yes you could only use one but since we both know $x$ and $y$ it is a
bit easier to use both than to only use one as you need to know $\theta$ at that point as well which may or may not be the case.}. So let's combine \autoref{eq:polar x} and \autoref{eq:polar y}
in \autoref{eq:polar lambda1}, transform it such that we end up with only $\lambda$ on one side and the rest on the other side and we end up with \autoref{eq:polar lambda3}.
\begin{equation}
\tan(\alpha) = \frac{\sin(\alpha)}{\cos(\alpha)}
\label{eq:tan}
\end{equation}
\begin{subequations}
\begin{equation}
\frac{x}{y} = \frac{a\sin(\lambda)}{a\cos(\lambda)} = \frac{\sin(\lambda)}{\cos(\lambda)}
\label{eq:polar lambda1}
\end{equation}
\begin{equation}
\frac{x}{y} = \tan(\lambda)
\label{eq:polar lambda2}
\end{equation}
\begin{equation}
\lambda = \tan^{-1}(\frac{x}{y})
\label{eq:polar lambda3}
\end{equation}
\end{subequations}
With this math we can fix a lot of stuff in the model. With this we can resample (mapping from sphere to plane) the pressure, density, temperarature and advection to the plane and ensure that
there are no more overflows and funky business. The implementation (code) for this will be done in a follow up stream, so stay tuned!

55
tex-docs/streams/TTNMETAF.tex
View File

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

116
tex-docs/topics/advection.tex Normal file
View File

@ -0,0 +1,116 @@
\section{Advection}
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.

123
tex-docs/topics/control_panel.tex Normal file
View File

@ -0,0 +1,123 @@
\section{Control Panel} \label{sec:cp}
Before we dive in an start modelling the planet, let us first set up a control panel that will influence how the model will behave and effectively decides what type of planet we model.
\subsection{The Beginning}
In the beginning there was nothing, and then there was "Hello World!" Or at least that is how many projects start. Why? you might ask, which is a perfectly valid question. In Computer Science,
"Hello World!" is very simple code that we use to test whether all the tools we need to get coding works. This checks whether the computer compiles the code and is able to execute it and whether
the code editor (IDE, Integrated Development Environment) starts the right processes to get the code compiled and executed. Oh right we were talking about CLAuDE, ahem.
Every project must have its beginning. And with CLAuDE I made the decision to start explaining the Control Panel first. This is to get you familiar with notation and to lay down some basics. To
do that we start with the fixed part of the Control Panel, the physical constants. Many things vary from planet to planet, how much radiation they receive from their star, how strong their
gravity is, how fast they spin around their axis and many many more. What does not change are the physical constants, well because they are constant. The Stefan-Boltzmann constant for instance
does not change. Whether you are on Earth, in space or on Jupiter, the value of the Stefan-Boltzmann constant will remain the same.
The Stefan-Boltzmann constant is denoted by $\sigma$ and has a value of $5.670373 \cdot 10^-8 \ (Wm^{-2}K^{-4})$ \cite{stefan-boltzmann}. The $\sigma$ is a greek letter called sigma. Greek
letters are often used in mathematics, as well as in physics or any other discipline that relies on maths (spoiler alert, quite a lot). Treat it like a normal letter in maths, representing a
number that you either do not know yet or is too long or cumbersome to write down every time. The Stefan-Boltzmann constant is denoted in scientific notation, a number followed by the order of
magnitude. It is denoted as a multiplication, because that is what you have to do to get the real number. An example: $4.3 \cdot 10^2 = 430$ and $4.3 \cdot 10^{-2} = 0.043$. The
letters behind the numbers are units, how we give meaning to the numbers. If I say that I am $1.67$ does not mean anything. Do I mean inches, centimeters, meters, miles? That is why we need units
as they give meaning to the number. they tell us whether the number is talking about speed, distance, time, energy and many other things. In this manual we will use SI units. Behind all the
letters you will find the following: [number]. This is a citation, a reference to an external source where you can check whether I can still read. If I pull a value out somewhere I will insert a
citation to show that I am not making these numbers up. This is what scientists use to back up their claims if they do not want to redo the work that others have done. I mean what is the point
of re-inventing the wheel if there is a tyre company next door? That is why scientists citate.
So with that out of the way, let us write down some constants. Why do I do this here? Because a lot of constants are used everywhere and I am too lazy to relicate them every time. If you see a
letter or symbol that is not explicitly explained, then it is most likely a constant that we discuss here in the control panel.
\subsection{Physical Constants}
As mentioned before, physical constants do not change based on where you are in the universe. Below you will find an overview of all the relevant constants together with their units. And a short
explanation where they are used or what they represent. To see them in action, consult the other sections of this manual, you will find them in equations all throughout this document.
\subsubsection{The Gas Constant}
The Gas constant, $R = 8.3144621$ ($J(mol)^{-1}K$) \cite{idealGas} is the constant used to relate the temperature of the gas to the pressure and the volume. One would expect this constant to be
different per gas, but under high enough temperatures and low enough pressure the gas constant is the same for all gases.
\subsubsection{The Specific Heat Capacity}
The specific heat capacity $c$ depicts how much energy is required to heat the object by one degree Kelvin per unit mass ($\frac{J}{Kg \cdot K}$) \cite{specificHeat}. This varies per material
and is usually indicated by a subscript. The specific heat capacity for water for instance is $c_w = 4190 JKg^{-1}K^{-1}$. Specific heat capacities also exist in the form of $Jg^{-1}K^{-1}$,
$Jmol^{-1}K^{-1}$ and $Jcm^{-3}K^{-1}$ which you can use in various circumstances, depending on what information you have.
\subsubsection{Mole}
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}. This is not a physical constant perse, but more like a unit ($mol$). Though it is still important enough to be added here for future reference. All other units are
way more intuitive and are assumed to be known.
\subsubsection{The Stefan-Boltzmann Constant}
The Stefan-Boltzmann constant, $\sigma = 5.670373 \cdot 10^-8 \ (Wm^{-2}K^{-4})$ \cite{stefan-boltzmann} is used in the Stefan-Boltzmann law (more on that in \autoref{sec:first thermolaw}).
\subsection{Planet Specific Variables}
The following set of variables vary per planet, that's why we call them variables since they vary. Makes sense right? We add them here as we will use them throughout the manual. The advantage
of that is quite significant. If you want to test things for a different planet, you only need to change the values in one place, instead of all places where you use it. If there is one thing
that we computer scientists hate is doing work, we like being lazy and defining things in one place means that we can be lazy if we need to change it. So we put in the extra work now, so we do
not have to do the extra work in the future. That's actually a quite accurate description of computer scientists, doing hard work so that they can be lazy in the future.
\subsubsection{The Passage of Time}
On Earth we have various indications of how much time has passed. While most of them remain the same throughout the universe, like seconds, minutes and hours, others vary throughout the universe,
like days, months and years. Here we specify how long the variable quantities of time are for the planet we want to consider as they are used in the code. This can be seen in
\autoref{alg:time constants}. Here a $\leftarrow$ indicates that we assign a value to the variable name before it, so that we can use the variable name in the code instead of the value, which
has the advantage I indicated before. // means that we start a comment, which is text that the code ignores and does not tell the cpu about. Not that the cpu would understand it, but that just
means less work for the computer. Yes computers are lazy too.
\begin{algorithm*}
\caption{Definition of how much time it takes for a day and a year on a planet and how much time on the planet passes before we start another calculation run}
\label{alg:time constants}
\SetKwComment{Comment}{//}{}
$day \leftarrow 60*60*24$ \Comment*[l]{Length of one day in seconds ($s$)}
$year \leftarrow 365*day$ \Comment*[l]{Length of one year in seconds ($s$)}
$\delta t \leftarrow 60 * 9$ \Comment*[l]{How much time is between each calculation run in seconds ($s$)}
\end{algorithm*}
\subsubsection{The Planet Passport}
Each planet is different, so why should they all have the same gravity? Oh wait, they don't. Just as they are not all the same size, tilted as much and their atmospheres differ. So here we define
all the relevant variables that are unique to a planet, or well not necessarily unique but you get the idea. This can all be found in \autoref{alg:planet constants}.
\begin{algorithm}
\caption{Defining the constants that are specific to a planet}
\label{alg:planet constants}
\SetKwComment{Comment}{//}{}
$g \leftarrow 9.81$ \Comment*[l]{Magnitude of gravity on the planet in $ms^{-2}$}
$\alpha \leftarrow -23.5$ \Comment*[l]{By how many degrees the planet is tilted with respect to the star's plane}
$top \leftarrow 50*10^3$ \Comment*[l]{How high the top of the atmosphere is with respect to the planet surface in meters ($m$)}
$ins \leftarrow 1370$ \Comment*[l]{Amount of energy from the star that reaches the planet per unit area ($Jm^{-2}$)}
$\epsilon \leftarrow 0.75$ \Comment*[l]{Absorbtivity of the atmosphere, fraction of how much of the total energy is absorbed (unitless)}
$r \leftarrow 6.4*10^6$ \Comment*[l]{The radius of the planet in meters ($m$)}
\end{algorithm}
\subsubsection{Model Specific Parameters}
These parameters cannot be found out in the wild, they only exist within our model. They control things like the size of a cell on the latitude longitude grid (more on that in later sections),
how much time the model gets to spin up. We need the model to spin up in order to avoid 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 it for a bit. You will
see that the whole system oscilates wildly as it adjusts to the sudden changes and eventually it will stabilise. We define the amount of time it needs to stabilise as the spin up time. All
definitions can be found in \autoref{alg:model constants}. What the $adv$ boolean does is enabling or disabling advection, a process described in <section> which does not work yet.
\begin{algorithm}
\caption{Defining the paramters that only apply to the model}
\label{alg:model constants}
\SetKwComment{Comment}{//}{}
$resolution \leftarrow 3$ \Comment*[l]{The amount of degrees on the latitude longitude grid that each cell has, with this setting each cell is 3 degrees latitude high and 3 degrees
longitude wide}
$nlevels \leftarrow 10$ \Comment*[l]{The amount of layers in the atmosphere}
$\delta t_s \leftarrow 60*137$ \Comment*[l]{The time between calculation rounds during the spin up period in seconds ($s$)}
$t_s \leftarrow 5*day$ \Comment*[l]{How long we let the planet spin up in seconds ($s$)}
$adv \leftarrow \texttt{FALSE}$ \Comment*[l]{Whether we want to enable advection or not}
$adv\_boun \leftarrow 8$ \Comment*[l]{How many cells away from the poles where we want to stop calculating the effects of advection}
$C_a \leftarrow 287$ \Comment*[l]{Heat capacity of the atmosphere in $JKg^{-1}K^{-1}$}
$C_p \leftarrow 1 \cdot 10^6$ \Comment*[l]{Heat capacity of the planet in $JKg^{-1}K^{-1}$}
$\delta y \leftarrow \frac{2\pi r}{nlat}$ \Comment*[l]{How far apart the gridpoints in the y direction are (degrees latitude)}
$\alpha_a \leftarrow 2 \cdot 10^{-5}$ \Comment*[l]{The diffusivity constant for the atmosphere}
$\alpha_p \leftarrow 1.5 \cdot 10^{-6}$ \Comment*[l]{The diffusivity constant for the planet surface}
$count \leftarrow 0$ \;
\For{$j \in [0, top]$}{
$heights[j] \leftarrow count$ \Comment*[l]{The height of a layer}
$count \leftarrow count + \frac{top}{nlevels}$ \;
}
\For{$i \in [0, nlat]$}{
$\delta x[i] \leftarrow \delta y\cos(lat[i]\frac{\pi}{180})$ \Comment*[l]{How far apart the gridpoints in the x direction are (degrees longitude)}
}
\For{$k \in [0, nlevels - 1]$}{
$\delta z[k] \leftarrow heights[k + 1] - heights$ \Comment*[l]{How far apart the gridpoints in the z direction are ($m$)}
}
\end{algorithm}

112
tex-docs/topics/master.tex Normal file
View File

@ -0,0 +1,112 @@
\section{The Master File}
The master file is the file that controls the model calculation. This file decides what calculations are used and what is done with the calculations (which is not the scope of this manual).
In other words, the master file combines all the calculations and theory from the previous sections and puts it all together to form a model. As mentioned earlier, this structure enables the
user to create their own version of the model. If one has their own calculations, or wants to use an older version of the calculations in this manual, then the user can define it themselves
and call it instead of the calls that we use. The model is meant to be customisable, which this structure enables.
\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:advectionfix} \;
}
}
}
\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.
\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{Non-uniform air density}
While air density on the surface is in general consistent, this does not hold if you move up through the atmosphere. The planet will pull air down due to gravity, which means that more air is at
the planet surface than at the top of the atmosphere. Hence the air density changes throughout the atmosphere and we need to account for that. This is done in \autoref{alg:density}. Because this
is used in radiation, velocity and advection, we initialise this in the master file. Though one could argue it could be part of the control panel, we opt not to include any code other than
variable declarations in the control panel for greater clarity. This also means that we give the user the option to have only one layer (by skipping implementing this algorithm). 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}
\subsection{Interpolating the Air Density}
In order to interpolate (see \autoref{sec:interpolation}) the air density, to have a better estimation at each grid cell, 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.

583
tex-docs/topics/radiation.tex Normal file
View File

@ -0,0 +1,583 @@
\section{Radiation} \label{sec:rad}
Radiation is energy waves, some waves are visible like light, others are invisible like radio signals. As is the basis for physics, energy cannot be created nor destroyed, only changed from one
form to another.
\subsection{The First Law of Thermodynamics and the Stefan-Boltzmann Equation} \label{sec:first thermolaw}
If energy goes into an object it must equal the outflowing energy plus the change of internal energy. Which is exactly what happens with the atmosphere. Radiation from the sun comes in, and
radiation from the atmosphere goes out. And along the way we heat the atmosphere and the planet which causes less radiation to be emitted than received. At least, that is the idea for Earth which
may not apply to all planets. Let one thing be clear, more radiation cannot be emitted than is inserted, unless the planet and atmosphere are cooling. Anyway, we assume that the planet is a black
body, i.e. it absorbs all radiation on all wavelengths. This is captured in Stefan-Boltzmann's law (\autoref{eq:stefan-boltzmann}) \cite{stefan-boltzmann}.
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 second per meters squared $Wm^{-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}
The energy difference between the energy 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 need to get the energy that we get ($J$) and not the energy per second ($W$), 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 function of the radiation part 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$), which is all done in \autoref{sec:cp}. 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. Do note that from here on most, if not all
functions need to be called \footnote{In case you are unfamiliar with calls, defining a function is defining how it works and calling a function is actually using it.} from another file which I
will call the master. The master file decides what parts of the model to use, what information it uses for plots and the like. We do it this way because we want to be able to switch out
calculations. Say that I find a more efficient way, or more detailed way, to calculate the temperature change. If everything was in one file, then I need to edit the source code of the project.
With the master file structure, I can just swap out the reference to the project's implementation with a reference to my own implementation. This makes the life of the user (in this case the
programmer who has another implementation) easier and makes changing calculations in the future easier as well. Also note that what we pass on as parameters \footnote{Parameters are variables
that a function can use but are defined elsewhere. The real values of the variables are passed on to the funciton in the call.} in \autoref{alg:stream1v1} are the
things that change during the execution of the model or that are calculated beforehand and not constants. $S$ for instance is not constant (well at this point it is but in \autoref{sec:daynight}
we change that) amd the current time is obviously not constant. All constants can be found in \autoref{sec:cp}.
\begin{algorithm}[hbt]
\SetAlgoLined
\SetKwInput{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{time $t$, amount of energy that hits the planet $S$}
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
\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}$ \;
}
}
\caption{The main function of the temperature calculations}
\label{alg:stream1v1}
\end{algorithm}
\subsection{Day/Night Cycle} \label{sec:daynight}
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. We need $-t$ as this ensures that the sun moves in the right direction, with $t$ the sun would move in the opposite direction in our model
than how it would move in real life. 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).
\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}.
\begin{algorithm}[hbt]
\SetAlgoLined
\SetKwInput{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{amount of energy that hits the planet $S$}
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
\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}$ \;
}
}
\caption{The main function of the temperature calculations}
\label{alg:stream1v2}
\end{algorithm}
\subsection{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.2. Do note that we change $C_p$ from a constant
to an array. We do this to allow adding in oceans or other terrain in the future. Same thing for the albedo, different terrain has different reflectiveness.
\begin{algorithm}[hbt]
$a \leftarrow 0.2$ \;
$C_p \leftarrow 10^6$ \;
\caption{Defining albedo}
\label{alg:albedo}
\end{algorithm}
Now that we have that defined, we need to adjust the main loop of the program (\autoref{alg:stream1v2}). 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
\SetKwInput{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{amount of energy that hits the planet $S$}
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
\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}$ \;
}
}
\caption{The main loop of the temperature calculations}
\label{alg:stream2v3}
\end{algorithm}
\subsection{Temperature with Varying Density}
The air density is not at all points exactly the same. This may be due to the winds blowing, or due to height changes in the terrain. We need to account for that, which is done in
\autoref{alg:temperature with density}.
\begin{algorithm}[hbt]
\SetAlgoLined
\SetKwInput{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{amount of energy that hits the planet $S$}
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
\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}$ \;
}
}
\caption{The main loop of the temperature calculations}
\label{alg:temperature with density}
\end{algorithm}
\subsection{Adding Layers}
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 calculation 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. 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$. 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}. $\delta z$
\begin{algorithm}[hbt]
\SetAlgoLined
\SetKwInput{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{amount of energy that hits the planet $S$}
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
\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]}$ \;
}
}
}
}
\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}
\subsection{Grey Radiation Scheme}
Inspired by the Isca project \cite{isca} and a paper describing the grey radiation scheme\cite{greyRad}.
A radiation scheme is a model for how energy is redistributed using light in a system. Such a model is a Grey radiation scheme if you split it into two parts, short and long wavelength radiation.
So you have two redistribution systems, one for short wavelength light and one for long wavelength light. Another assumption we make when using the Grey radiation scheme, is that the atmosphere
is transparent to short wavelength radiation, meaning it lets through light with short wavelengths. Additionally we use a two stream approximation, which means that we have a stream of radiation
going up, and another stream of radiation going down. Note that these two streams are both long wavelength radiation, because we said earlier we assume the atmosphere completely ignores short
wavelength radiation.
The two long wavelength radiation streams are described in \autoref{eq:upward radiation} and \autoref{eq:downward radiation} \cite{greyRad}. In those equations, the symbols are:
\begin{itemize}
\item $U$: Upward flux.
\item $D$: Downward flux.
\item $B$: The Stefan-Boltzmann equation (see \autoref{eq:stefan-boltzmann}).
\item $\tau$: Optical depth.
\end{itemize}
\begin{subequations}
\begin{equation}
\frac{dU}{d\tau} = U - B
\label{eq:upward radiation}
\end{equation}
\begin{equation}
\frac{dD}{d\tau} = B - D
\label{eq:downward radiation}
\end{equation}
\end{subequations}
With \autoref{eq:upward radiation} and \autoref{eq:downward radiation} written down, we can discuss how they work. These equations need a boundary condition to work, a starting point if you like.
For those equations the boundary conditions are that $U$ is at the surface equal to $B$ and that $D$ at the top of the atmosphere is equal to $0$. Meaning that in the beginning the top of the
atmosphere has no downward flux as there is no heat there, and that the bottom of the atmosphere has a lot of upward flux as most if not all of the heat is located there. Then after the spin up
time this should stabilise. We are interested in the change of the fluxes, so $dU$ and $dD$, to get those we need to multiply the right hand side by $d\tau$. Then we have the flow of radiation
that we need. However we cannot solely use these two equations to calculate the heat of a given layer. For that we need a few more components. These are described in \autoref{eq:heat layer}.
Here $Q_R$ is the amount of heat in a layer, $c_p$ is the specific heat capacity of dry air (our atmosphere), $\rho$ is the density of the air in that layer and $\delta z$ is the change in height.
$\delta U - D$ are the change in net radiation, meaning the amount of radiation that is left over after you transferred the upward and downward flux. See it as incoming and outgoing energy for a
given layer, the net change (either cooling down or heating up) is what remains after you have subtracted the incoming energy from the outgoing energy. While this explanation is not entirely true
(as flux is not entirely equivalent to energy), it explains the concept the best.
\begin{equation}
Q_R = \frac{1}{c_p\rho}\frac{\delta(U - D)}{\delta z}
\label{eq:heat layer}
\end{equation}
Now only one question remains: what is optical depth? Optical depth is the amount of work a photon has had to do to get to a certain point. This might sound really vague, but bear with me.
Optical depth describes how much stuff a certain photon has had to go through to get to a point. As you'd expect this is $0$ at the top of the atmosphere as space is a big vacuum so no stuff to
move through, so no work. Then the further the photon moves into the atmosphere, the more work the photon has had to do to get there. This is because it now needs to move through gases, like air,
water vapour and other gases. Hence the closer the photon gets to the surface of the planet, the larger the optical depth is because the photon has had to work more to get there. This phenomenon
is described in \autoref{eq:optical depth}. The symbols in the equation mean:
\begin{itemize}
\item $\tau_0$: Optical depth at the surface.
\item $p$: Atmospheric pressure ($Pa$).
\item $p_s$: Atmospheric pressure at the surface ($Pa$).
\item $f_l$: The linear optical depth parameter, with a value of 0.1.
\end{itemize}
\begin{equation}
\tau = \tau_0[f_l(\frac{p}{p_s}) + (1 - f_l)(\frac{p}{p_s})^4]
\label{eq:optical depth}
\end{equation}
As one can see, \autoref{eq:optical depth} has two parts, a linear part and a quatric part (to the power $4$). The quatric term approximates the structure of water vapour in the atmosphere, which
roughly scales with $\frac{1}{4}$ with respect to the height. The linear term is present to fix numerical behaviour because this is an approximation which will not be completely correct (that's
why it is an approximation) so we add this term to make it roughly right. The same thing holds for $f_l$ which can be manually tuned to fix weird numerical behaviour.
With these equations in our mind, let's get coding. First we add the pressure profile, the pressure of all atmospheric layers at a lat lon point. We need this to accurately represent the optical
depth per atmospheric layer. Then we need to use the pressure profile with regards to \autoref{eq:optical depth}. The resulting code can be found in \autoref{alg:optical depth}. This algorithm
replaces the temperature calculations we have done in \autoref{alg:temperature layer}, as this is basically a better version of the calculations done in that algorithm. $f_l$ has a value of $0.1$
and is located near all the other constants in the code, henceforth we will refer to this section in the code as the control panel, since most if not all of the constants can be tweaked here.
$\tau_0$ is a function that gives the surface optical depth for a given latitude. The corresponding equation can be found in \autoref{eq:optical depth surface} \cite{simon}. Translating this
into code is left as an exercise to the reader. $U[0]$ is the boundary condition discussed before (being the same as \autoref{eq:stefan-boltzmann}), just as $D[nlevels]$ is the boundary condition.
$S_z$ represents the call to \autoref{alg:gradient z}. \texttt{solar} represents the call to \autoref{alg:solar}.
\begin{algorithm}
\SetKwInput{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{amount of energy that hits the planet $S$}
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
\For{$lat \in [-nlat, nlat]$}{
\For{$lon \in [0, nlon]$}{
$pressureProfile \leftarrow p[lat, lon, :]$ \;
$\tau = \tau_0(lat)f_l\frac{pressureProfile}{pressureProfile[0]} + (1 - f_l)(\frac{pressureProfile}{pressureProfile[0]})^4)$ \;
$U[0] \leftarrow \sigma T_p[lat, lon]^4$ \;
\For{$level \in [1, nlevels]$}{
$U[level] \leftarrow U[level - 1] - \frac{(\tau[level] - \tau[level - 1])(\sigma \cdot (mean(T_a[:, :, level]))^4)}{1 + (\tau[level - 1] - \tau[level])}$ \;
}
$D[nlevels - 1] \leftarrow 0$ \;
\For{$level \in [nlevels - 1, 0]$}{
$D[level] \leftarrow D[level + 1] - \frac{(\tau[level + 1] - \tau[level])(\sigma \cdot (mean(T_a[:, :, level]))^4)}{1 + (\tau[level] - \tau[level + 1])}$ \;
}
\For{$level \in [0, nlevels]$}{
$Q[level] \leftarrow - \frac{S_z(U - D, 0, 0, level)}{10^3 \cdot densityProfile[level]}$ \;
}
$T_a[lat, lon, :] \leftarrow T_a[lat, lon, :] + Q$ \;
$S \leftarrow \texttt{solar}(I, lat, lon, t)$ \;
$T_p[lat, lon] \leftarrow T_p[lat, lon] \frac{\delta t((1 - a[lat, lon]) S + S_z(D, 0, 0, 0) - \sigma T_p[lat, lon]^4)}{C_p[lat ,lon]}$ \;
}
}
\caption{Adding in radiation}
\label{alg:optical depth}
\end{algorithm}
\begin{equation}
\tau_0 = 3.75 + \cos(lat \frac{\pi}{90})\frac{4.5}{2}
\label{eq:optical depth surface}
\end{equation}
\subsection{Adding In Some Ozone (Or Something Else That Approximates It)}
Adding in ozone in the stratosphere is hella complicated, so we leave that as an exercise to the reader as in true academic fashion. Just joking, if you want you can work on implementing ozone
however we opt not to because it is quite complicated. Instead we approximate it, which is decent enough for our purpose. We need to do it in \autoref{alg:optical depth} as we need to adjust the
$Q$. We add in a check to see if we are currently calculating the radiation in the stratosphere. If so we add some radiation extra to replicate the effect of ozone. As can be seen in
\autoref{alg:ozone}, where we only focus on the $Q$ part of \autoref{alg:optical depth}, we add in some extra radiation based on how high the current layer calculation is, which scales with the
height.
\begin{algorithm}
\For{$level \in [0, nlevels]$}{
$Q[level] \leftarrow - \frac{S_z(U - D, 0, 0, level)}{10^3 \cdot densityProfile[level]}$ \;
\uIf{$heights[level] > 20 \cdot 10^3$}{
$Q[level] \leftarrow Q[level] + \texttt{solar}(5, lat, lon, t) \frac{24 \cdot 60 \cdot 60(\frac{heights[level] - 20 \cdot 10^3}{10^3})^2}{30^2}$ \;
}
}
\caption{Replicating the effect of ozone}
\label{alg:ozone}
\end{algorithm}
\subsection{Tilting the Planet}
In order to model a planet that has seasons, like Earth, we need to tilt the planet. This has as effect that the sun is not always directly above the equator but is above a certain band around
the equator as the year moves on. This means that some hemispheres receive more/less sun based on what part of the year it is. Which corresponds to the various seasons we have on Earth. But in
order to do that, we have to change the \texttt{solar} function. The new version as shown in \autoref{alg:solar tilt} will replace \autoref{alg:solar}. Here $\alpha$ is the tilt in degrees.
\begin{algorithm}
\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, longitude and time combination.}
$sun\_lon \leftarrow -t \text{ mod } d$ \;
$sun\_lon \leftarrow sun\_lon \cdot \frac{360}{d}$ \;
$sun\_lat \leftarrow \alpha\cos(\frac{2t\pi}{year})$ \;
$S \leftarrow insolation\cos(\frac{\pi(lat - sun\_lat)}{180})$ \;
\uIf{$S < 0$}{
\Return{$0$} \;
} \uElse {
$lon\_diff \leftarrow lon - sun\_lon$ \;
$S \leftarrow S\cos(\frac{lon\_diff\pi}{180})$ \;
\uIf{$S < 0$}{
\uIf{$lat + sun\_lat > 90$ or $lat + sun\_lat < -90$}{
\Return{$insolation\cos(\frac{\pi(lat + sun\_lat)}{180})\cos(\frac{lon\_diff\pi}{180})$} \;
} \uElse {
\Return{$0$} \;
}
} \uElse {
\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 tilt}
\end{algorithm}
What the code in \autoref{alg:solar tilt} does boils down to calculating the latitude and longitude of the sun and checking whether the planet receives any energy. If not return $0$ immediately.
If so we check if the difference between the sun's longitude and the planet's longitude and calculate how much energy would hit the planet given that the sun is not straight above the equator.
We do this by multiplying the energy it would receive from the sun if it were above the equator $S$ by the cosine of the difference in longitudes, which represents the tilt. Then we check again
if the planet is receiving energy, if not we check if it happens around the poles. We do this because due to the tilt it can be the case that at certain points in the year the pole is in constant
sunlight, i.e. the sun does not go down. This creates a sort of overshoot which needs to be accounted for. If it does this then we add the latitudes of the sun and the planet together and use
that to calculate the energy that would hit that spot. If it is not the case that we are around the poles and we do not receive energy, then we return $0$. If it happens to be that we do receive
energy (so no negative values) then we return $S$.
\subsection{The Theory Behind the Planar Approximation}
It is time to deal with the pole situation. The north and south poles that is, not the lovely people over in Poland. We run into problems because the latitude longitude grid cells become to small
near the poles. Therefore, the magnitudes no longer fit into one cell and overflow into other cells which makes everything kind of funky. So we need to fix that, and we do that by a planar
approximation.
As said earlier, the grid cells on the latitude longitude grid get closer together the closer you get to the poles which poses problems. To fix this, we will be using a planar approximation of
the poles. What this means is that we will map the 3D grid near the poles onto a 2D plane parallel to the poles, as if we put a giant flat plane in the exact center of the poles and draw lines
from the grid directly upwards to the plane. For a visual representation, please consult the stream with timestamp 1:38:25 \cite{polarPlane}, which includes some explanation. In the streamm we
use $r$ to indicate the radius of the planet (which we assume is a sphere), $\theta$ for the longitude and $\lambda$ for the latitude. So we have spherical coordinates, which we need to transform
into $x$ and $y$ coordinates on the plane. We also need the distance between the center point (the point where the plane touches the planet which is the center of the pole) and the projected
point on the plane from the grid (the location on the plane where a line from the gird upwards to the plane hits it). This distance is denoted by $a$ (Simon chose this one, not me). We then get
the following equations as shown in \autoref{eq:polar distance}, \autoref{eq:polar x} and \autoref{eq:polar y}.
\begin{subequations}
\begin{equation}
a = r \cos(\theta)
\label{eq:polar distance}
\end{equation}
\begin{equation}
x = a \sin(\lambda)
\label{eq:polar x}
\end{equation}
\begin{equation}
y = a \cos(\lambda)
\label{eq:polar y}
\end{equation}
\end{subequations}
But what if we know $x$ and $y$ and want to know $\theta$ and $\lambda$? Pythagoras' Theorem then comes into play \cite{pythagoras}. We know that (due to Pythagoras) \autoref{eq:pythagoras} must
always be true. Then if we substitue $a$ by $\sqrt{x^2 + y^2}$ in \autoref{eq:polar distance} we get \autoref{eq:polar theta1}. Then we transform that equation such that we only have $\theta$ on
one side and the rest on the other side (since we want to know $\theta$) and we get \autoref{eq:polar theta3}.
\begin{equation}
x^2 + y^2 = a^2
\label{eq:pythagoras}
\end{equation}
\begin{subequations}
\begin{equation}
\sqrt{x^2 + y^2} = r\cos(\theta)
\label{eq:polar theta1}
\end{equation}
\begin{equation}
\frac{\sqrt{x^2 + y^2}}{r} = \cos(\theta)
\label{eq:polar theta2}
\end{equation}
\begin{equation}
\cos^{-1}(\frac{\sqrt{x^2 + y^2}}{r}) = \theta
\label{eq:polar theta3}
\end{equation}
\end{subequations}
For $\lambda$ we need another trigonometric function which is the tangent ($\tan$). The tangent is defined in \autoref{eq:tan}. If we then take a look at \autoref{eq:polar x} and
\autoref{eq:polar y}, we see that $\lambda$ is present in both equations. So we need to use both to get $\lambda$ \footnote{Yes you could only use one but since we both know $x$ and $y$ it is a
bit easier to use both than to only use one as you need to know $\theta$ at that point as well which may or may not be the case.}. So let's combine \autoref{eq:polar x} and \autoref{eq:polar y}
in \autoref{eq:polar lambda1}, transform it such that we end up with only $\lambda$ on one side and the rest on the other side and we end up with \autoref{eq:polar lambda3}.
\begin{equation}
\tan(\alpha) = \frac{\sin(\alpha)}{\cos(\alpha)}
\label{eq:tan}
\end{equation}
\begin{subequations}
\begin{equation}
\frac{x}{y} = \frac{a\sin(\lambda)}{a\cos(\lambda)} = \frac{\sin(\lambda)}{\cos(\lambda)}
\label{eq:polar lambda1}
\end{equation}
\begin{equation}
\frac{x}{y} = \tan(\lambda)
\label{eq:polar lambda2}
\end{equation}
\begin{equation}
\lambda = \tan^{-1}(\frac{x}{y})
\label{eq:polar lambda3}
\end{equation}
\end{subequations}
With this math we can fix a lot of stuff in the model. With this we can resample (mapping from sphere to plane) the pressure, density, temperarature and advection to the plane and ensure that
there are no more overflows and funky business. The implementation (code) for this will be done in a follow up stream, so stay tuned!

192
tex-docs/topics/util_funcs.tex Normal file
View File

@ -0,0 +1,192 @@
\section{Utility Functions}
With the control panel defined and explained, let us move over to some utility functions. Functions that can be used in all kinds of calculations, which we might need more often. In general it
concerns functions like calculating the gradient, the lacplacian or interpolation.
\subsection{Gradients}
Let us define the gradient in the $x, y$ and $z$ directions. The functions can be found in \autoref{alg:gradient x}, \autoref{alg:gradient y} and \autoref{alg:gradient z}. We use these functions
in various other algorithms as the gradient (also known as derivative) is often used in physics. It denotes the rate of change, how much something changes over time. Velocity for instance denotes
how far you move in a given time. Which is a rate of change, how much your distance to a given point changes over time.
In \autoref{alg:gradient z} $a.dimensions$ is the attribute that tells us how deeply nested the array $a$ is. If the result is $1$ we have just a normal array, if it is $2$ we have a double array
(an array at each index of the array) which is also called a matrix and if it is $3$ we have a triple array. We need this because we have a one-dimensional case, for when we do not use multiple
layers and a three-dimensional case for when we do use multiple layers. This distinction is needed to avoid errors being thrown when running the model with one or multiple layers.
This same concept can be seen in \autoref{alg:gradient x} and \autoref{alg:gradient y}, though here we check if $k$ is defined or \texttt{NULL}. We do this as sometimes we want to use this
function for matrices that does not have the third 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 third dimension in the gradient. We also use forward differencing
(calculating the gradient by taking the difference of the cell and the next/previous cell, multiplied by $2$ to keep it fair) in \autoref{alg:gradient y} as that gives better results for the
calculations we will do later on.
\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}
\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}
\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$}
\Output{Gradient in the $z$ direction}
\uIf{$a.dimensions == 1$}{
\uIf{$k == 0$}{
$grad \leftarrow \frac{a[k + 1] - a[k]}{\delta z[k]}$ \;
}\uElseIf{$k == nlevels - 1$}{
$grad \leftarrow \frac{a[k] - a[k - 1]}{\delta z[k]}$ \;
}\uElse{
$grad \leftarrow \frac{a[k + 1] - a[k - 1]}{2\delta z[k]}$ \;
}
} \uElse {
\uIf{$k == 0$}{
$grad \leftarrow \frac{a[i, j, k + 1] - a[i, j, k]}{\delta z[k]}$ \;
}\uElseIf{$k == nlevels - 1$}{
$grad \leftarrow \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]}{2\delta z[k]}$ \;
}
}
\Return $grad$ \;
\caption{Calculating the gradient in the $z$ direction}
\label{alg:gradient z}
\end{algorithm}
\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}
The code can be found in \autoref{alg:laplacian}. $\Delta_x$ and $\Delta_y$ in \autoref{alg:laplacian} represents the calls to \autoref{alg:gradient x} and \autoref{alg:gradient y}
respectively.
\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}
\end{algorithm}
\subsection{Divergence}
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 vectors $u, v$ and $w$ here already,
as we expect that we might use it in combination with the divergence operator more frequently. What those vectors are and represent we will discuss in \autoref{sec:momentum}.
\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, k] \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}
\end{algorithm}
\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}

157
tex-docs/topics/velocity.tex Normal file
View File

@ -0,0 +1,157 @@
\section{Air Velocity}
Did you ever feel the wind blow? Most probably. That's what we will be calculating here. How hard the wind will blow. This is noted as velocity, how fast something moves.
\subsection{Equation of State and the Incompressible Atmosphere}
The equation of state relates one or more variables in a dynamical system (like the atmosphere) 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 Momentum Equations} \label{sec:momentum}
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 movement, 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}
With the gradient functions defined in \autoref{alg:gradient x} and \autoref{alg:gradient y}, 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)$ \;
\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{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 replace it with \autoref{alg:coriolis}.
\begin{algorithm}
\SetAlgoLined
$\Omega \leftarrow 7.2921 \cdot 10^{-5}$ \;
\For{$lat \in [-nlat, nlat]$}{
$f[lat] \leftarrow 2\Omega \sin(lat \frac{\pi}{180})$ \;
}
\caption{Calculating the coriolis force}
\label{alg:coriolis}
\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}.
\subsection{Adding in Layers}
With adding in atmospheric layers 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 change in pressure. Therefore we need to store a copy of the pressure before we do any
calculations. This needs to be a copy due to aliasing \footnote{Aliasing is assigning a different name to a variable, while it remains the same variable. Take for instance that we declare a
variable $x$ and set it to be $4$. Then we say $y \leftarrow x$, which you might think is the same as saying they $y \leftarrow 4$ but behind the screen it is pointing to $x$. So if $x$ changes,
then so does $y$.}
\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)$ \;
\For{$lat \in [1, nlat - 1]$}{
\For{$lon \in [0, nlon]$}{
\For{$layer \in [0, nlevels]$}{
$u[lan, lon, layer] \leftarrow u[lat, lon, layer] + \delta t \frac{-u[lat, lon, layer]S_{xu} - v[lat, lon, layer]S_{yu} + f[lat]v[lat, lon, layer] - S_{px}}{\rho}$ \;
$v[lan, lon, layer] \leftarrow v[lat, lon, layer] + \delta t \frac{-u[lat, lon, layer]S_{xv} - v[lat, lon, layer]S_{yv} - f[lat]u[lat, lon, layer] - S_{py}}{\rho}$ \;
$w[lan, lon, layer] \leftarrow w[lat, lon, layer] - \frac{p[lat, lon, layer] - p_o[lat, lon, layer]}{\delta t\rho[lat, lon, layer]g}$ \;
}
}
$p_o \leftarrow copy(p)$ \;
}
\caption{Calculating the flow of the atmosphere (wind)}
\label{alg:velocity}
\end{algorithm}