mirror of https://github.com/Askill/claude.git
Finished the rewrite, only need to brush up the algorithms a bit
This commit is contained in:
TechWizzart
parent
e17e66b639
commit
47a20e14d6
|
|
@ -26,9 +26,16 @@
|
||||||
|
|
||||||
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
|
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,
|
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
|
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
|
||||||
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
|
Simon's Twitch which has all the streams without ad breaks or interruptions \cite{playlist}.
|
||||||
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}.
|
|
||||||
|
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
|
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
|
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
|
||||||
|
|
@ -43,24 +50,12 @@ This manual is for the toy model, which is as of now still in development. One i
|
||||||
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.
|
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
|
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
|
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
|
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).
|
knowing).
|
||||||
|
|
||||||
\input{streams/Stream5.tex}
|
|
||||||
|
|
||||||
\input{streams/Stream6.tex}
|
|
||||||
|
|
||||||
\input{streams/Stream7.tex}
|
|
||||||
|
|
||||||
\input{streams/Stream8.tex}
|
|
||||||
|
|
||||||
\input{streams/Stream9.tex}
|
|
||||||
|
|
||||||
\input{streams/Stream10.tex}
|
|
||||||
|
|
||||||
\input{topics/control_panel.tex}
|
\input{topics/control_panel.tex}
|
||||||
|
|
||||||
\input{topics/util_funcs.tex}
|
\input{topics/util_funcs.tex}
|
||||||
|
|
@ -74,7 +69,11 @@ knowing).
|
||||||
\input{topics/master.tex}
|
\input{topics/master.tex}
|
||||||
|
|
||||||
\newpage
|
\newpage
|
||||||
\input{streams/TTNMETAF.tex}
|
\input{appendices/TTNMETAF.tex}
|
||||||
|
|
||||||
|
\input{appendices/history.tex}
|
||||||
|
|
||||||
|
\input{appendices/vars.tex}
|
||||||
|
|
||||||
\newpage
|
\newpage
|
||||||
\bibliography{references}
|
\bibliography{references}
|
||||||
|
|
|
||||||
|
|
@ -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).
|
||||||
|
|
@ -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}
|
||||||
|
|
@ -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}
|
|
||||||
}
|
|
||||||
|
|
@ -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).
|
|
||||||
|
|
@ -1,144 +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}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
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}
|
|
||||||
|
|
||||||
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$.
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
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}
|
|
||||||
|
|
@ -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}
|
|
||||||
|
|
@ -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.
|
|
||||||
|
|
@ -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.
|
|
||||||
|
|
@ -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!
|
|
||||||
|
|
@ -111,4 +111,28 @@ boundary.
|
||||||
}
|
}
|
||||||
\caption{The main loop for calculating the effects of advection}
|
\caption{The main loop for calculating the effects of advection}
|
||||||
\label{alg:advectionfix}
|
\label{alg:advectionfix}
|
||||||
\end{algorithm}
|
\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}
|
||||||
|
$\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}
|
||||||
|
|
||||||
|
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.
|
||||||
|
|
@ -43,7 +43,7 @@ mass units ($u$)\cite{mole}. This is not a physical constant perse, but more lik
|
||||||
way more intuitive and are assumed to be known.
|
way more intuitive and are assumed to be known.
|
||||||
|
|
||||||
\subsubsection{The Stefan-Boltzmann Constant}
|
\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 %insert reference here).
|
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}
|
\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
|
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
|
||||||
|
|
@ -80,7 +80,7 @@ all the relevant variables that are unique to a planet, or well not necessarily
|
||||||
$top \leftarrow 50*10^3$ \Comment*[l]{How high the top of the atmosphere is with respect to the planet surface in meters ($m$)}
|
$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}$)}
|
$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)}
|
$\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$)}
|
$r \leftarrow 6.4*10^6$ \Comment*[l]{The radius of the planet in meters ($m$)}
|
||||||
\end{algorithm}
|
\end{algorithm}
|
||||||
|
|
||||||
\subsubsection{Model Specific Parameters}
|
\subsubsection{Model Specific Parameters}
|
||||||
|
|
@ -101,4 +101,21 @@ definitions can be found in \autoref{alg:model constants}. What the $adv$ boolea
|
||||||
$t_s \leftarrow 5*day$ \Comment*[l]{How long we let the planet spin up 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 \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}
|
$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)}
|
||||||
|
|
||||||
|
$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}
|
\end{algorithm}
|
||||||
|
|
@ -27,8 +27,7 @@ calculations in \autoref{alg:temperature with density} we would calculate the ve
|
||||||
|
|
||||||
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
|
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
|
\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
|
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.
|
||||||
recommend leaving it on \texttt{FALSE} until it is fixed in \autoref{sec:advectionfix}.
|
|
||||||
|
|
||||||
\begin{algorithm}
|
\begin{algorithm}
|
||||||
\eIf{$t < 7day$}{
|
\eIf{$t < 7day$}{
|
||||||
|
|
@ -40,4 +39,57 @@ recommend leaving it on \texttt{FALSE} until it is fixed in \autoref{sec:advecti
|
||||||
}
|
}
|
||||||
\caption{The spin-up block dynamically enabling or disabling flow simulation}
|
\caption{The spin-up block dynamically enabling or disabling flow simulation}
|
||||||
\label{alg:spinup}
|
\label{alg:spinup}
|
||||||
\end{algorithm}
|
\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.
|
||||||
|
|
@ -1,8 +1,8 @@
|
||||||
\section{Radiation}
|
\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
|
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.
|
form to another.
|
||||||
|
|
||||||
\subsection{The First Law of Thermodynamics and the Stefan-Boltzmann Equation}
|
\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
|
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
|
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
|
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
|
||||||
|
|
@ -214,12 +214,12 @@ to an array. We do this to allow adding in oceans or other terrain in the future
|
||||||
\begin{algorithm}[hbt]
|
\begin{algorithm}[hbt]
|
||||||
$a \leftarrow 0.2$ \;
|
$a \leftarrow 0.2$ \;
|
||||||
|
|
||||||
$C_p \leftarrow 10^7$ \;
|
$C_p \leftarrow 10^6$ \;
|
||||||
\caption{Defining the oceans}
|
\caption{Defining albedo}
|
||||||
\label{alg:albedo}
|
\label{alg:albedo}
|
||||||
\end{algorithm}
|
\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
|
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
|
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.
|
reflected instead of absorbed, where we need the amount that is absorbed which is exactly equal to $1$ minus the amount that is reflected.
|
||||||
|
|
||||||
|
|
@ -240,6 +240,9 @@ reflected instead of absorbed, where we need the amount that is absorbed which i
|
||||||
\end{algorithm}
|
\end{algorithm}
|
||||||
|
|
||||||
\subsection{Temperature with Varying Density}
|
\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]
|
\begin{algorithm}[hbt]
|
||||||
\SetAlgoLined
|
\SetAlgoLined
|
||||||
\SetKwInput{Input}{Input}
|
\SetKwInput{Input}{Input}
|
||||||
|
|
@ -272,4 +275,320 @@ being picked \cite{uniformdist}.
|
||||||
}
|
}
|
||||||
\caption{Varying the albedo of the planet}
|
\caption{Varying the albedo of the planet}
|
||||||
\label{alg:albedo variance}
|
\label{alg:albedo variance}
|
||||||
\end{algorithm}
|
\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
|
||||||
|
\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}
|
||||||
|
|
||||||
|
\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}
|
||||||
|
\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!
|
||||||
|
|
@ -41,145 +41,6 @@ important enough to account for yet, especially considering the current complexi
|
||||||
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
|
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}.
|
$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}
|
|
||||||
\textbf{NOTE:} This whole subsection is obsolete. We have replaced these calculations with \autoref{sec:momentum}. The folloing subsection is left in for historical value, and maybe for a simpler
|
|
||||||
calculation if you want your own model to do less heavy calculations. This is where the previously mentioned master file strucutre comes in. You can create a new file with the following
|
|
||||||
calculations and replace the call that you would make to \autoref{sec:momentum} with a call to the algorithm listed in this subsection. Your choice, though the model Simon has made opted to use
|
|
||||||
the more complicated calculations. So here are the original calculations and if you want an up to date overview of the calculations please have a look at \autoref{sec:momentum}.
|
|
||||||
|
|
||||||
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).
|
|
||||||
|
|
||||||
\subsection{The Momentum Equations} \label{sec:momentum}
|
\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
|
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
|
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
|
||||||
|
|
@ -266,4 +127,35 @@ Another change introduced is in the coriolis parameter. Up until now it has been
|
||||||
|
|
||||||
\subsection{Adding Friction}
|
\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
|
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}.
|
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)$ \;
|
||||||
|
\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}
|
||||||
Loading…
Reference in New Issue