claude/tex-docs/streams/Stream10.tex

\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).