\documentclass[10pt]{article}
\usepackage{a4wide}
\textwidth 18.5cm
\oddsidemargin -1.5cm
\textheight 25.5cm
\topmargin -2cm
\pagestyle{empty}


\begin{document}

\def\abs#1{\left| #1 \right|}
\def\localdir{.}


\title{Numerical modeling, Tutorial 2 : One dimensional advection of a tracer
\author{
Fr\'ed\'eric Hourdin,
hourdin@lmd.ens.fr
}\date{\today}}

\def\grad{\vec{\mbox{grad}}\ }
\def\afaire#1{{\bf A FAIRE : #1}}
\maketitle


Let us consider a trace species whose concentration $q(x,t)$, in kg/kg for instance, depends on dimension $x$ and time only, over a domain  $[0,L]$ in $x$ and $[0,T]$ in time.
Let us assume that this trace species is advected by a non divergent
 $U=$cste wind field, with a constant mass of air per unit length ($\rho$ in kg/m).
You may choose  $L=1000$~km and $U=10$~m~s$^{-1}$ to fix ideas or the legend of your plots, but those parameters are not important for the exercise.

\begin{enumerate}


\item  We assume that at initial time, the tracer is distributed 
as a Gaussian function around $x_0=L/5.$ , i.~e. $q(x,t=0)= exp( - [(x-x_0)/\lambda]^2 )$ with $\lambda=L/15.$
What will be the tracer concentration at a time $t$ after.

\item Given a discretization of 100 points in the domain $[0,L]$, write a fortran code which prints in column the values of  $x_i$ and
$q_{i,0}=q(x_i,t=0)$ for the 100 grid points.
Plot these values on the screen.
One can start from the program available at
\begin{verbatim}
http://www.lmd.jussieu.fr/~hourdin/COURS/M2/TD/TD2maquette.f90
\end{verbatim}
or on
\begin{verbatim}
~/MON/TD/TD2maquette.f90
\end{verbatim}

In this simple case, one can easily compute the  exact concentration
at time $t$.
Print the values at $t=0.6*L/U$ in a second file.
Superimpose the two curves on a graph.

Please parameterize the number of points in $x$, by declaring 
for instance an array for $x_i$ such as:

\begin{verbatim}
implicit none
integer, parameter :: imax=100
integer i

real, dimension(imax) :: xxx
do i=1,imax
   xxx(i)=...
\end{verbatim}


\item Give the 3 advection schemes based on 3 finite difference
estimates of the space derivative of the tracer concentration.
We will call upstream, centered and downstream schemes the
schemes centered respectively on a left, centered and right space
derivative (for $U>0$).
Based on a Taylor expansion in the vicinity of point $x_i$,
compare the accuracy of the various discretization schemes for the tracer concentration space derivative when $\delta x$ goes to 0. Which scheme is the most accurate ?

\item Write a time loop so as to compute the tracer advection
with an upstream scheme over 350 time steps, with $\delta t= 0.2 \times \delta x / U$.
Be careful of well separating in the code: 1) the declarations of variables, then  2) the initialization (attribution of a value) and
3) the computation itself.
Initialization must be done in particular before the time loop.\\
$[(x-x_0)/y]^2\rightarrow$ ((x-x0)/y))**2 in Fortran.
Also note that you should not have more than one time loop in your code, with as many lines as you want in x inside.


In order to record the results, one can write:

\begin{verbatim}
character*4 file
do it=1,nt
....
    if (mod(it,10)==1) then
       file="q..."
       write(file(2:4),'(i3.3)') it  ! writing, with 3 charcers, an intefer
                                     ! in position 2 to 4 of a
                                     ! 4-character string.
       open (10,file=file,form='formatted')
       do i=1,imax
          write(10,*) xxx(i),qqq(i)
       enddo
       close(10)
    endif
\end{verbatim}
which will create a file each 10 time step, called for instance
'q091' at time step 91.
Then curves can be plotted on the screen with the following command:
\begin{verbatim}
xmgrace q*
\end{verbatim}
You can also send the image in a file directly with the line command
\begin{verbatim}
xmgrace -hardcopy -hdevice JPEG -printfile out.jpg q*
\end{verbatim}
('q*' means that the command is applied for all the files starting by letter 'q' in your current directory, ie here : q001, q011, q021 ...).
But you can use any other mean to plot the curves.

The coding can be facilitated by introducing two "arrays" to store
the tracer concentrations $q$ at time 
$n$ and $n+1$, for instance 'qn(:)' and 'qnp1(:)'.

For the left boundary condition, one can assume $q=0$ for $x<0$.

\item  Adapt the program to use a centered scheme in $x$.
Compare the results of those two schemes.
Comment in terms of accuracy, by looking in particular at the first time steps, and with the Taylor expansion done above.

Comment on the more or less physical characteristics of the scheme.


\def\abs#1{\left| #1 \right|}
 
\def\localdir{.}



\def\xmax{x_{\mbox{max}}}
\def\zmax{z_{\mbox{max}}}
\def\mod{\mbox{mod}}

\item
Compute the advection with both schemes assuming that at initial time, the
tracer is the combination of two sin modes, and of a "numerical oscillation" at the grid scale :
 $q_{i,1}= \sin(4\pi x_i/\xmax) + \sin(17\pi x_i/\xmax) + OSC(i)/2$
where $OSC(i)=1$ for even and $-1$ for odd values of $i$, or, en fortran
\begin{verbatim}
pi=2.*asin(1.)
xxx(i)=i*deltax
qinitial(i)=sin(4*pi*xxx(i)/xmax)+sin(17*pi*xxx(i)/xmax)+0.5*(-1.+2.*mod(i,2))
\end{verbatim}
One may assume that flux entering on the left side or outgoing on the right
side of the domain are nuls, or alternatively code periodic boundary conditions.


\item Show that the upstream scheme can be written as the sum of the centered
scheme plus a term which can be interpreted as a finite difference formulation
of a diffusion operator
\begin{equation}\label{eq:diff}
\frac{\partial q}{\partial t} = K \frac{{\partial}^2 q}{\partial x^2}
\end{equation}
Give the value of K.
Apply the diffusion alone to the initial tracer field of the previous question.


\item Which scales dissipate first ? Comment with respect to the
time evolution of the amplitude $A(t)$ of a sin mode $A(t)\sin(k x)$ solution of Eq~\ref{eq:diff}.

\item Rerun the computation of dissipation over 10 time steps, saving results at each individual time step, but with  $\beta$=0.3 and 0.6, where $\beta = K \delta t / (\delta x)^2$.
Analyse the results obtained considering a solution of the discretized
scheme of the form $q_{i,n}=A_n OSC(i)$, where $n$ is the time step  (a link can be established with the results of tutorial \#1).


\end{enumerate}


Write a report on Tutorial 1 and 2.

Answer :\\
for TD1 : "Analysis" (1 page 2 figures max) \\
for TD2 : questions 3 and 5 to 9 (2 pages max including figures).

Once your repport is completed, save it as a PDF file. Name it as
\verb+TD12_name1_name2.pdf+, \verb+name1+  being your name (very important for me not to mix the various files). Then send the pdf file by e-mail at \verb+hourdin@lmd.jussieu.fr+.
{\bf put also your name and first name at the top of the document}.



\end{document}
