Subversion Repositories pvinversion.ecmwf

\documentclass[12pt]{article}

\renewcommand{\topfraction}{0.93}

\renewcommand{\bottomfraction}{0.93}

\renewcommand{\textfraction}{.07}

\input BoxedEPS

\SetRokickiEPSFSpecial

\ForceOn

\SetEPSFDirectory{M:/Weiterbildung/Allgemeine_Informatik/Diplomarbeit/fig/}

\HideDisplacementBoxes

\textwidth16.0cm

\textheight26.0cm

\oddsidemargin0.5cm

\topmargin-1.5cm

\headsep0cm

\topskip0cm

\renewcommand{\baselinestretch}{1.2}

\begin{document}

\pagenumbering{Roman}

% ----------------------------------------------------------------------------------------------

% Title page

% ----------------------------------------------------------------------------------------------

\thispagestyle{empty}

\vspace{0.7cm}

\begin{center}

{\Large\bf Numerical Piecewise Potential Vorticity Inversion \\[2mm]

A user guide for real-case experiments}

\vspace{1cm}

\begin{center}

\begin{minipage}[t]{10cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{Fig05_pv_320K.eps}

\end{minipage}

\end{center}

\vspace{1cm}

Michael Sprenger, Ph.\,D.$^{1}$ \\

\vspace{0.3cm}

{\sl Institute for Atmospheric and Climate Science, ETH Z\"urich, Switzerland}\\

\vspace{1.8cm}

\today \\

{Diploma Thesis for the Postgraduate Course in Computer Science, FHSS Schweiz,\\

under the supervision of Lars Kruse, Ph.\,D.}

\end{center}

\vspace{2cm}

\noindent

{\sl $^1$Corresponding author and his address:}\\[1mm]

Michael Sprenger\\

Institute for Climate and Atmospheric Science, {\it IAC}ETH \\

ETH Zurich \\

CH-8092 Zurich, Switzerland \\

e-mail: michael.sprenger@env.ethz.ch

\newpage

\vspace{20cm}

\noindent

{\bf Title figure:} Ertel's potential vorticity at 320\,K. 

% ----------------------------------------------------------------------------------------------

% Abstract

% ----------------------------------------------------------------------------------------------

\newpage

\begin{center}

{\large \bf Abstract}  

\end{center}

\noindent

Potential vorticity (PV) is a key quantity in atmospheric dynamics. Its value is based upon several

points: (a) PV is conserved in purely adiabatic flows, i.e. if radiative and condensational heating and other

diabatic processes can be neglected; (b) under suitable balance conditions the specification of PV

in the interior of a domain and with boundary values for potential temperature completely determines

the flow field, in particular horizontal wind and temperature can be derived; (c) in many respects

the atmosphere can be partitioned into several distinct PV features, which then interact and allow

the dynamics to be interpreted as the interaction of these PV elements. Points (b) and (c) are known

as the invertibility and partitioning principle, respectively.\\

\noindent

In this thesis, a program package is presented which allows to isolate PV elements and then to

study their impact on the atmospheric flow field and on the temperature distribution. Technically,

this is done with a so-called PV inversion, which in turn comprises several different steps, as for 

example transformation into suitable co-ordinate systems and numerical solution of a poisson equation. With PV inversion, two major problems arise: Firstly, a suitable balance condition must

be formulated, secondly the inversion problem for Ertel's PV is nonlinear, and therefore poses a significant challenge. The program presented in this thesis is based upon the so-called

quasi-geostrophic balance condition. More specifically, the linear quasi-geostrophic PV equation is

solved. Since quasi-geostrophic and Ertel's PV do not exactly coincide, an iterative technique

is adopted to approach the nonlinear Ertel-PV inversion by means of successively quasi-geostrophic

inversions.\\ 

\noindent

The aim of this work is to present an in-depth discussion of all steps which are needed for a PV inversion.

Special focus was given to a clear and user-friendly program package, where all steps are 

well documented and hence are attractive for further development. In this respect, the complete

 inversion is controlled by one single Linux Shell script.

The new user needs only to adjust some few paths in this script. The

main parameters for the inversion itself are specified in a separate parameter file.  

\newpage

\begin{center}

{\large \bf Zusammenfassung}  

\end{center}

\noindent

Die potentielle Vortizit\"at (PV) ist eine Kerngr\"osse der dynamischen Meteorologie. Ihre Bedeutung basiert auf mehreren Eigenschaften: (a) In einer adiabatischen Str\"omung, dh. wenn Strahlung,

die Freisetzung von latenter W\"arme und andere diabatische Prozesse vernachl\"assigt werden 

k\"onnen, ist PV eine lagrange'sche Erhaltungsgr\"osse; (b) unter geeigneten Balancebedingungen bestimmt die Verteilung der PV im Innern und der potentiellen Temperatur am Rand eines Gebietes

vollst\"andig den Zustand der Atmosph\"are, insbesondere die horizontalen Windfelder und das

Temperaturfeld; (c) h\"aufig l\"asst sich die Atmosph\"are in mehrere, klar voneinander getrennte PV-Elemente unterteilen, die dann die Dynamik der Atmosph\"are durch ihre Wechselwirkung bestimmen. 

Die Eigenschaften (b) und (c) sind in der Literatur bekannt als das Invertibilit\"atsprinzip und 

Partitionierungsprinzip.\\

\noindent

In dieser Arbeit wird ein Programmpaket vorgestellt, mit dessen Hilfe sich einzelne PV-Elemente

isolieren und ihr Einfluss auf die atmosph\"arischen Windfelder und Temperaturen studieren 

lassen. Programmtechnisch besteht diese sogenannte PV-Inversion aus mehreren Schritten. So muss

zum Beispiel eine Koordinatentransformation durchgef\"uhrt werden und eine Poisson-Gleichung

numerisch gel\"ost werden. Damit die PV-Inversion sinnvoll ist, m\"ussen zwei Probleme gel\"ost

werden. Zun\"achst muss eine geeignete Balancebedingung vorgegegeben sein. Ausserdem handelt es

sich bei der Inversion um einen nichtlinearen Prozess, der numerisch einige Herausforderungen

stellt. In dem Programmpaket dieser Arbeit wird die sogenannte quasi-geostrophische Balance

verwendet, dh. der mathematische/numerische Inversionsprozess besteht in der L\"osung der

linearen quasi-geostrophischen PV-Gleichung. Das Problem der Nichtlinearit\"at wird dadurch gel\"ost, dass

die Berechnung des genannten linearen Inversionsproblems mehrmals wiederholt wird. Iterativ

ergibt sich somit eine L\"osung des nichtlinearen Problems.\\

\noindent

Ziel dieser Arbeit ist eine detailierte Darstellung aller Schritte, die bei einer PV-Inversion

n\"otig sind. Besonderes Augenmerk wurde bei der Entwicklung des Programmpakets auf die klare

Strukturierung und Anwenderfreundlichkeit gelegt. Dadurch kann das Paket einerseits als Black-Box

verwendet werden, zugleich erlaubt es erfahrenen Benutzern eine leichte Einarbeitung in die 

Programme und damit die M\"oglichkeit zu deren Erweiterung. Die ganze Inversion wird von einem einzigen

Linux Shell Skript kontrolliert, in dem lediglich einige wenige Pfade angepasst werden m\"ussen. Die 

meisten Parameter, welche das gestellte Problem der PV-Inversion beschreiben, sind in einer

separaten Parameterdatei eingetragen.

% ----------------------------------------------------------------------------------------------

% Inhalt

% ----------------------------------------------------------------------------------------------

\newpage

\section*{Contents}

\vspace{0.5cm}

\renewcommand{\contentsname}{}

\contentsline{section}{\hspace{0.5cm} Abstract}{III}

\contentsline{section}{\hspace{0.5cm} Zusammenfassung}{IV}

\contentsline{section}{\hspace{0.5cm} Abbreviations}{VII}

\contentsline{section}{\hspace{0.5cm} List of Figures}{VIII}

\vspace{-1.5cm}

\tableofcontents

% ----------------------------------------------------------------------------------------------

% Abkürzungen

% ----------------------------------------------------------------------------------------------

\newpage

\section*{Abbreviations}

\vspace{0.5cm}

\begin{tabular}{rll}

{\bf [1]} & {\bf DWD}           & Deutscher Wetterdienst (German weather service)\\[2mm]

{\bf [2]} & {\bf ECMWF}         & European Centre for Medium Range Weather Forecasts\\[2mm]

{\bf [3]} & {\bf ERA-40}        & Re-analysis of the ECMWF for the periods 1958-2001\\[2mm]

{\bf [4]} & {\bf ERA-15}        & Re-analysis of the ECMWF for the periods 1979-1993\\[2mm]

{\bf [5]} & {\bf ETH}           & Eidgen\"ossische Technische Hochschule\\[2mm]

{\bf [6]} & {\bf f}            & Coriolis parameter (in $s^{-1}$)\\[2mm]

{\bf [7]} & {\bf IACETH}        & Institute for Atmospheric and Climate Science at ETH Zurich\\[2mm]

{\bf [8]} & {\bf METEOSAT}      & European geostationary weather satellite\\[2mm]

{\bf [9]} & {\bf MeteoSwiss}    & Swiss weather service\\[2mm]

{\bf [10]} & {\bf NaN }          & Not a Number (missing data flag)\\[2mm]

{\bf [11]} & {\bf netcdf     }   & Compact and random-access file format\\[2mm]

{\bf [12]} & {\bf NWP     }      & Numerical weather prediction\\[2mm]

{\bf [13]} & {\bf $N^2$    }     & Squared Brunt-Vais\"al\"a frequency (in $s^-2$)\\[2mm]

{\bf [14]} & {\bf PV}            & Potential vorticity (if not otherwise stated Ertel's potential vorticity)\\[2mm]

{\bf [15]} & {\bf pvu}            & Potential vorticity unit\\[2mm]

{\bf [16]} & {\bf PDE}           & Partial differential equation\\[2mm]

{\bf [17]} & {\bf SOR}           & Successive over-relaxation\\[2mm]

{\bf [18]} & {\bf T}             & Temperature (in $^\circ$C)\\[2mm]

{\bf [19]} & {\bf $T_v$}         & Virtual temperature (in $^\circ$C)\\[2mm]

{\bf [20]} & {\bf U}             & Zonal velocity component (in west/east direction)\\[2mm]

{\bf [21]} & {\bf V}             & Meridional velocity component (in south/north direction)\\[2mm]

{\bf [22]} & {\bf $\omega$}      & Vertical velocity component (in Pa/s)\\[2mm]

{\bf [23]} & {\bf q}             & Specific humidity (in kg/kg)\\[2mm]

{\bf [24]} & {\bf $\theta$}      & Potential temperature (in K)\\[2mm]

{\bf [25]} & {\bf $\rho$}        & Air density (in kg/m$^3$)\\[2mm]

{\bf [26]} & {\bf R$_d$}         & Gas constant for dry air (in kg/m$^3$)\\[2mm]

{\bf [27]} & {\bf $\epsilon$}    & Ratio of gas constants for dry air and water vapour\\[2mm]

{\bf [28]} & {\bf g}             & Earth's gravity (in kg m/s$^2)$\\

\end{tabular}

\newpage

\listoffigures

% ----------------------------------------------------------------------------------------------

% Introduction and Motivation

% ----------------------------------------------------------------------------------------------

\newpage

\clearpage

\pagenumbering{arabic}

\section{Introduction and Motivation}

\begin{figure}[t]

\begin{center}

\begin{minipage}[t]{16cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{sat_wv.ps}

\end{minipage}

\\

\caption{\it Water vapour satellite imagery (in gray shading) and isolines of Ertel's potential vorticity (PV)

at 350\,hPa and for 3 May 1998, 12\,UTC. Dark regions indicate a dry upper troposphere, whereas bright and white regions are

indicative for a lot of moisture in the upper troposphere. Note that high-PV regions are coincident with 

dry regions. }

\label{satwv}

\end{center}

\end{figure}

There are different ways how the state of the atmosphere can be described. Traditionally this is done by specification of temperature, horizontal and vertical wind components and of either pressure, if geometrical height is used as the vertical co-ordinate, or geopotential height, if pressure is used as the vertical co-ordinate. This traditional approach has its main advantage in its simplicity. On the other hand, theoretical physics has often experienced the case that new insight can be gained if abstract, but physically more fundamental quantities were introduced. Consider for example the introduction of action and its primary physical parameter, the Planck constant. Similarly, geophysical fluid dynamics and hence atmospheric dynamics made some transitions away from the simple traditional meteorological fields toward more abstract ones.\\

\noindent

A long existing concept of engineering fluid dynamics is the concept of vorticity, which mathematically is defined as the curl of the wind vector field. Vorticity is treated in nearly every text book on fluid dynamics, for instance in Acheson (1990). Vorticity is a scalar quantity, but nevertheless, in a divergent-free and two-dimensional flow its specification is sufficient to deduce the complete velocity field. We could call this the invertibility principle for vorticity in such a non-divergent and two-dimensional flow. Moreover, an interesting conservation law hold for vorticity under these assumptions: Following the fluid motion, vorticity is conserved, i.e. its lagrangian derivative with respect to time is zero. Probably, in its most elegant way this conservation principle is expressed in Kelvin's circulation theorem (Acheson, 1990). These concepts of barotropic flow and of conservation of vorticity were indeed the basis for the first numerical weather predictions by Charney in 1950 (for an interesting history of numerical weather prediction, consult either Lynch, 2006, or Nebeker, 1995).\\ 

\noindent

Of course, the atmosphere cannot be treated in an exact way as a barotropic fluid. It constitutes a

stratified fluid where density and pressure decreases with increasing height. Therefore, a suitable 

generalisation of the barotropic-vorticity equation for the real atmosphere is by far not trivial, but would be highly desirable. 

Potential vorticity goes into this direction. The fundamental work by Rossby (1940) and Ertel (1942)

showed that potential vorticity is conserved in adiabatic (no latent heat release, no radiative heating

or cooling) and frictionless flow. The generalisation of the "invertibility principle for vorticity"

in two-dimensional flow was first stated in a rough way by Kleinschmidt (1950). He was able to attribute

some low-level flow features to an upper-level PV anomaly, in his words to a "Zyklonenk\"orper". Figure\,1

is a nice illustration of an upper-level PV anomaly. It shows some isolines of Ertel's PV on 350\,hPa, 

which form an elongated and narrow filament of high PV extending from Scandinavia south to the north-western

edge of Spain. Remarkably, this anomaly of high PV is also discernible as a dark band in the water vapour 

imagery of the geostationary METEOSAT weather satellite. Hence, the high PV band is associated with a very dry upper troposphere. \\

\begin{figure}[t]

\begin{center}

\begin{minipage}[t]{16cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{Fig11_PVInv_tropo1.ps}

\end{minipage}

\\

\caption{\it Section across the circularly symmetric structure induced by an isolated, circularly symmetric, cyclonic potential vorticity anomaly near the model tropopause across which the Rossby-Ertel potential vorticity has a strong discontinuity, by a factor of 6. The bold line is the tropopause, the horizontal

thin lines are potential temperature, plotted at 5\,K intervals, and the circular/elliptical lines denote tangential velocity, at 3\,m/s intervals. The structure is typical for middle latitudes; the Coriolis parameter is $10^{-4}s^{-1}$ (as at geographical latitude $43.3^\circ$N). The domain shown has a radius of 2500\,km (taken from McIntyre,

1997).}

\label{inversion}

\end{center}

\end{figure}

\noindent

A very influential set of equations was introduced by Charney by means of a scale analysis of the dynamical equations (1948). These so called quasi-geostrophic equations are well suitable to describe synoptic and planetary-scale

processes, whilst neglecting smaller-scale features. Although these equations nowadays are considered no

longer of sufficient accuracy for numerical weather prediction, they nevertheless still are of great importance in

theoretical dynamic meteorology due to their simplicity and elegance (Holton, 1992).

Along with this set of equations came a new version of potential vorticity. This quasi-geostrophic

potential vorticity is linearly "linked" to the flow and temperature field, and -particularly interesting for

this study- an invertibility principle can be formulated. Indeed, the linear relationship between quasi-geostrophic potential vorticity and the flow field makes it very interesting since sophisticated

numerical techniques exist for the solution of this kind of problem. An example of such an inversion is shown in Fig.\,\ref{inversion} where

a cyclonic potential-vorticity anomaly near a model tropopause is shown. The potential vorticity is high

in the stratosphere and low in the troposphere. Therefore, the downward excursion of the tropopause is 

associated locally with anomalously high potential vorticity. This anomaly is marked in the figure with

stippling. The "horizontal" lines correspond to the potential temperature, the "circular/elliptical" contours to the tangential velocity. Note that the isolines of potential temperature (the so-called isentropes) are

pulled upwards below the PV anomaly, and pulled downward within the PV anomaly. An interesting aspect of

the upward pulling of the isentropes is the associated reduction in atmospheric stability. Due to the reduced

vertical gradient of  potential temperature, the atmosphere is more prone to convective instability. In fact,

in a recent study Martius et al. (2006) showed that many heavy precipitation events in the south of the Alps are 

associated with such upper-level distortions of the tropopause, and hence with a PV anomaly. Moreover, note

that the upper-level PV anomaly is not only associated with a deformation of the isentropes, but is

also linked with a wind field. The induced wind field in this idealised setting reaches 21\,m/s at its maximum and

is cyclonically (anti-clockwise) oriented. The wind field is strongest at the tropopause level, but is even discernible at

the surface far below the anomaly. In this respect, the PV field exhibits a far-field effect, and this in turn

is the basic idea behind the invertibility principle. If this principle is taken together with the partitioning 

principle, its explaining power becomes particularly attractive. This latter principles states that the atmospheric state can

be expressed as the interaction of isolated PV elements. This immediately invokes a research strategy, which is: to

enter, modify or remove some of these PV elements, and to see how the flow evolves in the thus

changed state. A very influential review article on PV and PV thinking was presented by Hoskins et al. (1985).

Here the key elements of PV thinking are discussed and applied to many atmospheric dynamical problems. \\

\begin{figure}[t]

\begin{center}

\begin{minipage}[t]{16cm}

  \vspace{-2cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{streamer.eps}

\end{minipage}

\\[-2cm]

\caption{\it Stratospheric PV streamer and positions with stratosphere-troposphere mass exchange. The section is

for 1 January 1986, 06\,UTC and on the 320\,K isentrope. It is taken from the ERA-15 data set of the ECMWF. In color Ertel's potential vorticity is shown in pvu. Mass exchange from the stratosphere to the troposphere is

marked with stars, transport in the opposite direction by open circles (taken from Sprenger et al., 2007)}

\label{streamer}

\end{center}

\end{figure}

\noindent

Many studies relate to this PV thinking perspective. For instance, Davis and Emanuel (1991) looked at the 

potential vorticity diagnostics of cyclogenesis; Fehlmann and Davies (1997) investigated the impact of PV 

structures on precipitation events over Europe. In addition to this methodology, PV has gained a lot of

interest in recent years as a diagnostic tool. Indeed, from a dynamicist's point of view it is highly

attractive to use potential vorticity as the defining component of the extra-tropical tropopause (Stohl et al. 2003). Typically,

this value is set to 2\,pvu, and every air parcel with larger PV is treated as stratospheric and

every air parcel with lower PV as tropospheric. This definition allows to discuss the dynamics of the extra-tropical tropopause, as it would not be feasible with the more common thermal (or lapse-rate) definition of the tropopause. In Fig.\,\ref{streamer} a prominent feature of the extra-tropical tropopause is shown: a stratospheric

PV streamer. It corresponds to a pronounced excursion of stratospheric (high-PV) air towards the south. These

features are quite common in the extra-tropics and are important for the mass exchange between the stratosphere

and the troposphere, amongst other impacts (for example they are often related to cyclogenesis). The positions where stratospheric mass (and chemical constituents, as for example ozone) crosses the tropopause are marked in the figure by stars. They are predominantly found on the upstream

(western) side of the streamer. Crossing in the opposite direction, i.e. from the troposphere to the stratosphere,

on the other hand occurs on the downstream (eastern) side of the streamer. Interesting questions emerge from the

inspection of a such a figure: If the intensity, or the southward excursion of the PV streamer was reduced, would there still be a significant mass flux across the tropopause? These, and similar questions are at the heart of

the PV inversion, which is introduced in this study.\\

\noindent

The study is organised in the following way: In chapter 2 the mathematical problem of PV inversion is 

formulated and some essential aspects of the numerical algorithm are discussed. Chapter~3 introduces some

key aspects for the re-structuring of the program package.

Then, in chapter 4 a detailed

example of PV inversion is discussed. This part is intended to be a practical user guide for PV inversion,

and therefore rather detailed technical aspects are presented.  Chapter 5 follows with some helpful 

diagnostic tools which allow to quality and impact of the PV inversion. Finally, chapter 6 concludes 

with some general remarks and wishes regarding the PV inversion tool.

\newpage

\section{Mathematical Framework and Numerical Aspects}

In this study we use the quasi-geostrophic approximation for the real-case

inversion problem. In this limit the relative (quasi-geostrophic) potential vorticity

$q$ is given by:\\[-2mm]

\[

q = \frac{\partial^2\psi}{\partial x^2} + \frac{\partial^2\psi}{\partial y^2} +

\frac{f^2}{\rho_0} \cdot \frac{\partial}{\partial z} ( \frac{\rho_0}{N_0^2}

\cdot \frac{\partial \psi}{\partial z})

\]

\vspace{0.3cm}

\noindent

with the boundary values of potential temperature at the upper and lower lid and 

of zonal and meridional wind components at the lateral boundaries:\\[-2mm]

\[

g \cdot \frac{\theta^*}{\theta_0} = f \cdot \frac{\partial \psi}{\partial z}

\quad \quad

u = -\frac{\partial \psi}{\partial y}

\quad \quad

v  =  -\frac{\partial \psi}{\partial x} 

\]

\vspace{0.3cm}

\noindent

Here $N_0^2$, $\rho_0$ and $\theta_0$  denote the squared Brunt-V\"as\"ala frequency, the density and the

potential temperature of a reference state depending only on the vertical co-ordinate $z$. $\psi$ is the 

streamfunction from which the horizontal wind components can be derived. Finally, $f$ is the Coriolis parameter

which measures the Earth's rotation rate (typically $10^{-4}s^{-1}$ in the mid-latitudes and 0 at the equator),

and $g$ is the Earth's gravity.\\

\noindent

The above equations constitute a so-called von Neumann boundary problem for the elliptic

differential equation relating the potential vorticity $q$ to the streamfunction $\psi$. Hence,

the elliptic partial differential equation (PDE) can be transformed into the general form:\\[-2mm]

\[

\frac{\partial}{\partial x}(\alpha \cdot \frac{\partial \psi}{\partial x}) +

\frac{\partial}{\partial y}(\beta  \cdot \frac{\partial \psi}{\partial y}) +

\frac{\partial}{\partial z}(\gamma \cdot \frac{\partial \psi}{\partial z}) +

\sigma = 0

\]

\vspace{0.3cm}

\noindent

where the coefficients $\alpha$, $\beta$, $\gamma$ and $\sigma$ are easily expressed in terms of

the density $\rho_0$, the Coriolis parameter $f$ and the squared Brunt-Vais\"al\"a frequency $N^2$:\\[-2mm]

\[

\alpha = \rho_0(z) \quad \quad 

\beta  = \rho_0(z) \quad \quad 

\gamma = \frac{f^2(y) \cdot \rho_0(z)}{N_0^2(z)} \quad \quad 

\sigma = - \rho_0(z) \cdot q(x,y,z)

\]

\vspace{0.3cm}

\noindent

The next step is to discretise the above PDE by means of finite differences. We assume that all

fields (quasi-geostrophic PV $q$, streamfunction $\psi$,...) are defined on a three-dimensional grid, whose

grid points can be addressed by the three indices $i$, $j$, and $k$ in the x-, y- and z-direction (see Fig.\,ref{grid}). With this

convention, the discretised PDE can be written as:\\[-2mm]

\[

\Delta_x(A \cdot \Delta_x \psi) + 

\Delta_y(B \cdot \Delta_y \psi) +

\Delta_z(C \cdot \Delta_z \psi) +

S = 0

\]

\vspace{0.3cm}

\noindent

Here the operators $\Delta_x(A \cdot \Delta_x\psi)$, $\Delta_y(B \cdot \Delta_y\psi)$ and

$\Delta_z(C \cdot \Delta_z\psi)$ at the grid position $i,j,k$ are defined by:

\begin{eqnarray*}

\Delta_x(A \cdot \Delta_x\psi)_{i,j,k} & = &

A_{i+1/2,j,k} \cdot (\psi_{i+1,j,k}-\psi_{i,j,k}) -

A_{i-1/2,j,k} \cdot (\psi_{i,j,k}-\psi_{i-1,j,k}) \\[0.2cm]

\Delta_y(B \cdot \Delta_y\psi)_{i,j,k} & = &

B_{i,j+1/2,k} \cdot (\psi_{i,j+1,k}-\psi_{i,j,k}) -

B_{i,j-1/2,k} \cdot (\psi_{i,j,k}-\psi_{i,j-1,k}) \\[0.2cm]

\Delta_z(C \cdot \Delta_z\psi)_{i,j,k} & = &

C_{i,j,k+1/2} \cdot (\psi_{i,j,k+1}-\psi_{i,j,k}) -

C_{i,j,k-1/2} \cdot (\psi_{i,j,k}-\psi_{i,j,k-1})

\end{eqnarray*}

\vspace{0.3cm}

\noindent

Some simple algebra yields the following expressions for $A$, $B$ and $C$:

\begin{eqnarray*}

A_{i,j,k} & = & \frac{\Delta y \cdot \Delta z}{\Delta x} \cdot \alpha_{i,j,k}\\[0.2cm]

B_{i,j,k} & = & \frac{\Delta y \cdot \Delta z}{\Delta x} \cdot \beta_{i,j,k}\\[0.2cm]

C_{i,j,k} & = & \frac{\Delta y \cdot \Delta z}{\Delta x} \cdot \gamma_{i,j,k}

\end{eqnarray*}

\begin{figure}[t]

\begin{center}

\begin{minipage}[t]{16cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{numericalgrid.eps}

\end{minipage}

\\

\caption{\it The numerical grid for the PV inversion in an xz cross-section. The horizontal grid spacing is $\Delta x$

in the x-direction (and correspondingly $\Delta y$ in the y-direction), the vertical grid spacing $\Delta z$. Note that in the vertical, a staggered grid at half-levels is also needed (taken from Fehlmann, 1997).}

\label{grid}

\end{center}

\end{figure}

\vspace{0.3cm}

\noindent

Finally, the additive operator $S$ is expressed as:

\[

S_{i,j,k}  = \Delta x \cdot \Delta y \cdot \Delta z \cdot \sigma_{i,j,k}

\]

\vspace{0.3cm}

\noindent

It turns out to be advantageous to have the coefficients $\alpha$, $\beta$, $\gamma$ and $\sigma$

not only on the 3d grid where the quasi-geostrophic PV and the streamfunction are defined, but to

have them also on intermediate points. These leads to the grid structure which is shown in Fig.\,\ref{grid}.

On the left side the indices for the $\psi$-grid is shown, i.e. for the grid where the 

quasi-geostrophic PV and the streamfunction is defined. The right side gives the indices for the

coefficients $\alpha$, $\beta$, $\gamma$ and $\sigma$, i.e. for the intermediate layers.\\

\noindent

Let $\rho(k) = \rho_{k/2}$ and $N^2(k)=N^2_{k/2}$, hence expressing the vertical grid with intermediate layers. With these definitions the coefficients $A$, $B$ and $C$ for the inversion problem are readily obtained (note that these coefficients are dependent only on the vertical index z):

\begin{eqnarray*}

A_k & = & \frac{\Delta y \cdot \Delta z}{\Delta x} \cdot \rho(2k) \quad (k=0,...,nz)\\[0.2cm]

B_k & = & \frac{\Delta y \cdot \Delta z}{\Delta x} \cdot \rho(2k) \quad (k=0,...,nz)\\[0.2cm]

C_k & = & \frac{\Delta y \cdot \Delta z}{\Delta x} \cdot \frac{\rho(k) \cdot f^2}{N^2(k)} \quad (k=0,...2\cdot nz)

\end{eqnarray*}

\vspace{0.3cm}

\noindent

The additive operator $S$ depends on all three grid indices. In its discretised form it is:

\[

S_{i,j,k} = - \Delta x \cdot \Delta y \cdot \Delta z \cdot \rho(2k) \cdot q_{i,j,k} 

\quad (i=0,...,nx, j=0,...,ny, k=0,...,nz)

\]

\vspace{0.3cm}

\noindent

It is convenient to additionally define the coefficients $C_{-1}=C_0$ and $C_{2nz+1}=C_{2nz}$. With 

all these definitions the fully discretised quasi-geostrophic PV equation can be written as:

\begin{eqnarray*}

S_{i,j,k} & = &

       \quad \quad (\psi_{i-1,j,k}-2\psi_{i,j,k}+\psi_{i+1,j,k}) \cdot A_k       \\[0.2cm]

&  & + \quad (\psi_{i,j-1,k}-2\psi_{i,j,k}+\psi_{i,j+1,k}) \cdot B_k       \\[0.2cm]

&  & + \quad (\psi_{i,j,k-1} - \psi_{i,j,k})               \cdot C_{2k-1}  \\[0.2cm]

&  & + \quad (\psi_{i,j,k+1} - \psi_{i,j,k})               \cdot C_{2k+1} 

\end{eqnarray*}

\vspace{0.3cm}

\noindent

with the indices ranges as follows: $i=0, ... nx$ and $k=0, ... nz$. In addition to this discretised 

PV equation, additional equations must be specified for terms like $\psi_{-1,j,k}$. These discretised

von Neumann boundary conditions are:

\begin{eqnarray*}

A_k \cdot (\psi_{0,j,k} - \psi_{-1,j,k}) & = & {\Delta y} \cdot {\Delta z} \cdot \rho(2k) \cdot v_a(j,k)\\[0.2cm]

A_k \cdot (\psi_{nx+1,j,k} - \psi_{nx,j,k}) & = & {\Delta y} \cdot {\Delta z} \cdot \rho(2k) \cdot v_b(j,k)\\[0.2cm]

B_k \cdot (\psi_{i,0,k} - \psi_{i,-1,k}) & = & -{\Delta x} \cdot {\Delta z} \cdot \rho(2k) \cdot u_a(i,k)\\[0.2cm]

B_k \cdot (\psi_{i,ny+1,k} - \psi_{i,ny,k}) & = & -{\Delta x} \cdot {\Delta z} \cdot \rho(2k) \cdot u_b(i,k)

\end{eqnarray*}

\vspace{0.3cm}

\noindent

where $u_a,u_b$ denote the velocity components in x direction at the left (western) and right (eastern) lateral boundary, and

correspondingly  $v_a,v_b$ denote the velocity components in y direction at the front (southern) and back (northern) lateral boundary. Additional boundary values are specified at the lower and upper lid of the domain:

\begin{eqnarray*}

C_{-1} \cdot (\psi_{i,j,0} - \psi_{i,j,-1}) & = & {\Delta x} \cdot {\Delta z} \cdot 

\frac{f \cdot g \cdot \rho(0) \cdot \theta_{bot}(i,j)}{N^2(0) \cdot \theta_0(0)}\\[0.2cm]

C_{2 nz + 1} \cdot (\psi_{i,j,nz+1} - \psi_{i,j,nz}) & = & {\Delta x} \cdot {\Delta z} \cdot 

\frac{f \cdot g \cdot \rho(2 nz) \cdot \theta_{top}(i,j)}{N^2(2 nz) \cdot \theta_0(2 nz)}

\end{eqnarray*}

\vspace{0.3cm}

\noindent

Here, $\theta_{bot}$ and $\theta_{top}$ denote the boundary conditions for potential temperature at the 

lower and upper lid of the domain.\\

\noindent

This is a system of linear equations which can be expressed as\\[-2mm]

\[

B\psi = b

\]

\vspace{0.3cm}

\noindent

where $B$ is a $m \times m$ matrix with in total $m=(nx+1) \cdot (ny+1) \cdot (nz+1)$ elements (not to be 

confused with the above coefficients $B_k$). This linear

system has a solution if the following condition is fulfilled:

\[

\sum_{i,j,k} b_{i,j,k} = 0

\]

\vspace{0.3cm}

\noindent

This is exactly the discrete version of the compatibility condition in section~5.1. The derivation of the 

condition starts with the observation that the null space, i.e. the kernel of the system is at least 

one-dimensional because the non-trivial vector $\psi_{i,j,k}=1 \forall i,j,k$ is element of this kernel. Physically, this expresses the fact that the streamfunction is determined only up to an additive constant. Since the kernel of the linear system is at least one-dimensional, the image of the operator $B$ is

at most (m-1) dimensional. Moreover, the operator $B$ is normal, and therefore the image is orthogonal to the

kernel. Because $b$ is in the image of the operator $B$ and $\psi_{i,j,k}=1$ is in the kernel, the two vectors

are orthogonal. This leads immediately to the necessary consistency condition $\sum_{i,j,k} b_{i,j,k} = 0$.

Note that this expresses a complicated relationship which must be fulfilled by the interior PV distribution and the boundary values of potential temperature and horizontal wind components, because $b$ includes all these

forcing terms.\\

\noindent

The consistency condition is not met if vanishing boundary values are specified. This automatically

leads to an inconsistency in the numerical solution of the equation. In practice, the fulfillment

of the consistency condition can be enforced if a potential temperature "correction" is added at the

lower and upper lid, this correction being uniform and of opposite sign on the two boundaries. For 

reasonable PV and temperature distributions this additive temperature shift remains smaller than about 2\,K.\\

\noindent  

There exist several techniques how to solve linear systems of equations. Here we adopt the successive

over-relaxation (SOR) method: Let $A$ be a linear operator which is represented by an $m \times m$ matrix, and

let $\omega$ be a real number (the relaxation parameter) such that $|1-\omega (1+A)|< 1$. Then a solution

of the equation can be obtained iteratively, starting with an arbitrary first guess $\psi_i^{(0)}$. The iterations

\[

\psi_i^{(n+1)} = \omega \cdot (b_i - \sum_{j=1}^{i-1} A_{i,j} \cdot \psi_j^{(n+1)}

                                  - \sum_{j=i}^{m}   A_{i,j} \cdot \psi_j^{(n)})

                                  + (1-\omega) \cdot \psi_i^{(n)}

\]

\vspace{0.3cm}

\noindent

converge toward the solution of the system $A\psi + \psi =b$. If we choose $A=B-1$, the iterations converge toward a solution of the quasi-geostrophic PV equation. \\

\noindent

Note that the above outlined algorithm allows to overwrite the variable $\psi_i^{(n)}$ with the updated

variable $\psi_i^{(n+1)}$, and therefore needs a minimum of computational memory. In the Fortran program 

{\em inv\-cart.f} the number of iterations and the SOR parameter $\omega$ are specified. The are set to 500 iterations and 1.81, respectively.

% ----------------------------------------------------------------------------------------------

% Program philosophy 

% ----------------------------------------------------------------------------------------------

\newpage

\section{Re-Structuring of the Code}

In this section some key concepts of the program package are described. In fact, an earlier version of the package already existed (developed originally by Rene Fehlmann and later modified by Sebastien Dirren) and the present work is based upon this pre-existing package. So, the question arises what added value this re-development and re-structuring of the code is associated with. The basic 

idea is to provide a code which fulfills the following three requirements:

\begin{center}

\begin{minipage}[t]{10cm}

\begin{itemize}

\item[a)] {\bf Transparency and Modularity}\\[-3mm]

\item[b)] {\bf Universality and Model Independence}\\[-3mm]

\item[c)] {\bf User friendliness}

\end{itemize}

\end{minipage}

\end{center}

\vspace{0.3cm}

\noindent

In the following parts, a detailed discussion is intended to illustrate how this requirements were missing in the existing code

and how they could be incorporated into the new version. Before doing so, it is worthwhile to consider why the original code was lacking these requirements. A short review of the relevant literature

illustrates that the situation is actually quite common. Any software which is really used has to be adapted to the changing demands. It must be maintained. Unfortunately this maintenance often leads

to a degeneration of the code. The changes introduced into the code then make necessary some further

adaptions, and so on. Finally, in the end the code becomes so "distorted" and "chaotic" that often

a complete re-coding is easier to be done than a re-structuring of the existing version. The

remedy against this code degeneration is a continuous re-structuring. Hence, if changes in the code

are obligatory due to new user demands, it is important not only to introduce the needed changes

short-sightedly. The overall structure of the program should be kept in mind, and if possible

long-term perspectives in code maintenance should be allowed, although such a long-term perspective

momentarily leads to an additional effort. \\

\noindent

Computer science, in particular software engineering, has defined the term "re-structuring" (or

"re-factoring" for object-oriented programming languages) for the process how code should be

maintained in order to avoid severe degeneration. Fowler defines the term re-factoring in the

following way (taken from URL {\em www.refactoring.com}):\\

\begin{center}

\parbox{14cm}{"Refactoring is a disciplined technique for restructuring an existing body of code, altering its internal structure without changing its external behavior. Its heart is a series of small behavior preserving transformations. Each transformation (called a 'refactoring') does little, but a sequence of transformations can produce a significant restructuring. Since each refactoring is small, it's less likely to go wrong. The system is also kept fully working after each small refactoring, reducing the chances that a system can get seriously broken during the restructuring."}

\end{center}

\vspace{0.6cm}

\noindent

Hence, the key aspect of re-structuring is to continuously maintain the code, not only in

its performance, but also in its readability, documentation and adaptabilty to needed 

changes.\\

\vspace{0.3cm}

\noindent

 Scientific programming in particular is very prone to code degeneration. This is certainly due

 to the fact that cutting-edge science is, by "definition", investigating what is not known. 

 Therefore, if for instance a computer program is written to numerically simulate a physical 

 process, it might well turn out during the development of the code that complete new aspects

 must be included. This naturally leads to chaotically structured code. In the present case, i.e.

 for the PV inversion program, this certainly explains to a large part why the code became

 so un-organised. Several works have been done with the program, several researchers included

 into the code what they needed for their specific task, without considering that other users might  have no use of their changes. The results was a program package which, in principle, was still

 working, but the knowledge of how to apply it went lost.\\

 \noindent

  A first rough inspection immediately made

 clear that major changes were necessary to make the code available to a broader community again. 

 Would it have been possible to only re-organise the code and write a new user guide? In a first

 try, this was in fact considered. But then the advantages of a more or less complete

 re-coding turned out to be more suitable. Hence, the original software package joined the

 sad fate of so many degenerated codes: a complete re-coding. On the other hand, such a re-coding

 must be seen as a great chance. It allows to improve the code in such a way, as it would never be

 possible if only "slight" changes at the existing code were made.\\ 

 \noindent

 In the following three 

 subsections, some key aspects of this re-coding will be presented. They are by no means

 exhausting, but should give an impression of the new "philosophy". Hopefully, they also motivate

 future users of the program package to "successfully" maintain it.

\subsection{Transparency and Modularity}

What makes a computer program readable? Probably one of the most important points is "modularity".

The problem to be solved numerically can and should be split into several distinct steps. For instance,

for the PV inversion a classic three step splitting is appropriate: pre-processing, PV inversion, and post-processing. Moreover, each of these three primary steps can further be split into several secondary sub-steps. 

In the existing code, the splitting of the problem into distinct sub-problems was not clearly discernible. Indeed, the main Fortran program included not only the numerical inversion of the quasi-geostrophic PV equation, but also some of the

preparatory steps and some of the post-processing steps. \\

\noindent

A major improvement was gained by a very strict separation of the distinct primary steps. So, pre-processing, inversion and post-processing are done by three completely separated program packages. This strict 

separation is supported by the fact that three separated directories are used: There is one directory where pre-processing is done, on directory where inversion is done, and finally one directory where post-processing is done. A flow of data between the three directories is allowed only at well-defined

steps within the whole process. At the end of pre-processing the relevant files, and only those, are moved to the inversion (or run) directory. Similarly, at the end of the inversion, the relevant files

are moved from the run directory to the post-processing directory.\\ 

\noindent

Additional improvement resulted from a very clear separation of the sub-processes in the three main-processes (pre-processing, inversion, post-processing). In fact, it was not tried to incorporate

all preparatory steps into one large Fortran program, but instead each well-defined task (as for example transformation into a new co-ordinate system) constitutes a separate Fortran program. This 

modularity is very similar to the "philosophy" of the Linux operating system: Keep flexibility and 

clarity by offering not one program which handles everything (and thereby becomes a "monstrosity"), but offer many flexible and simple tools which the user only has to combine in order to perform complex tasks.\\

\noindent 

What else except for modularity can be done to improve the readability of computer programs? It is obvious, and probably the most neglected aspect of good computer code generation: in-code documentation. Every 

small sub-section of a computer program should be documented. It should be possible to gain 

insight into an algorithm only by looking at the in-code documentation. It is definitely not

sufficient only to document for every subroutine what its aim is and what its interface is, although

quite often even this most basic principle is not fulfilled. Hence, in the re-coding of the

PV inversion focus was led to good documentation: Where are files opened, where variables initialised, what is a subroutine call meaning? \\

\noindent

Finally a remark concerning the used programming language. A complete re-coding of a software package might also be

a chance to switch to a "better" programming language. In the present case, the original code was written in 

Fortran 77, and this language was also used for the re-coding. Fortran, and in particular its older version

Fortran 77, is by far not a modern language. Comparing its vocabulary and its structural elements to languages like C++ or Java, its power is very limited. Moreover, Fortran includes jump statements like {\em Goto} which could easily make

any code un-readable. So why do so many research codes still rely on Fortran? In fact, most numerical

weather prediction models are written in Fortran. There are two reasons why the PV inversion was written in Fortran: Firstly, the small vocabulary and the intuitive naming of the inherent Fortran commands

makes it quite easy to read a code. Many researches never had a profound introduction into programming. Nevertheless they should be able to implement algorithms for their daily work. This is easily done with Fortran, and certainly much easier done than using object-oriented concepts offered by C++ and other modern languages. The PV inversion should be available

to a research community which has excellent algorithmic thinking, but only a "weak" background in computer languages. Fortran is an optimal choice in this respect. The second reason why the re-coding was done in Fortran is its high performance. In fact, PV inversion is very resource demanding, and compilers must create fast-running codes in order

the PV inversion tool to be helpful.\\ 

\noindent

Regarding the "unholy" {\em Goto} statements, it can be immediately said that good programming style is not restricted to languages which do not offer a {\em Goto} statement. A well-structured and carefully developed Fortran program is

certainly preferable to a "bad" algorithm in a "good" computer language, "good" meaning that this language offers many

controlling structures. Probably, a key concept of well-written computer code is a good balance between the "power" of a (Fortran)

subroutine and its length (expressed in the number of code lines). Two extreme examples might illustrate this point: If no subroutines at all are used, i.e. if the program consists only in one single main program, the algorithm might get

"lost" in two many "technical" details. On the other hand, if every single and minor step is a separate subroutine, the

code lacks clarity due to the large number of subroutines. In the present re-coding of the PV inversion, focus was given

to a good balance of functionality (power) of a subroutine and its length. 

\subsection{Universality and Model Independence}

Universality was another aspect which was important in re-coding the PV inversion. In the present

example universality might be better named "model independence". Model in this respect refers either to different numerical weather prediction (NWP) models or to idealised experiments. In the former case, 

the variety of NWP models is very large. The PV inversion tool was originally written for the Europa

model (EM) of the German weather service (DWD). Later it was applied to the higher-resolution 

model (HRM and CHRM) which was run by the DWD and the Swiss weather service (MeteoSwiss). In recent years, the 

non-hydrostatic local model (LM) developed by the DWD replaced the HRM in regional

weather forecasting. So, an adaption of the existing PV inversion code would have been necessary. 

Moreover, PV inversion is particularly attractive for global-scale data sets. So, the method and the

programs were adapted to conform with the global model of the European Centre for Medium Range Forecasts (ECMWF). To make things worse, idealised experiments should also be performed. What actually

makes the situation so difficult is that each model has its own horizontal grid structure and its own vertical grid structure. For instance, ECMWF uses a terrain-following co-ordinate system which in higher levels becomes a pressure-level co-ordinate system, quite in contrast what is used by the regional LM model (a variant of geometrical height). \\

\noindent

Model independence cannot completely avoided. It is a fact which has to be adopted and cleverly 

dealt with. In the re-coding, model aspects were removed in the first two preparatory steps (out of 

eight different steps, see section~4.4) and the way back to the model is done in the last post-processing steps.

Essentially the preparatory steps interpolate the needed meteorological fields onto a stack of height levels and

 transforms the fields into a local quasi-cartesian co-ordinate system. After completing these steps, the other preparatory steps and the PV inversion itself is completely independent of the model from which the input data is retrieved. As a particular example consider the PV inversion of a structure over the North Pole. In the existing code such an inversion would not have been feasible due to the convergence of the

 longitude circles at the North Pole. The new code, on the other hand, elegantly circumvents this

 problem by introducing a new quasi-cartesian co-ordinate system centered over the North Pole. In this

 respect the North Pole is not different from any other region on the globe.\\

\noindent

With respect to idealised experiments, the following strategy was adopted. This kind of experiments 

is generally so different from what is needed for real-case studies that a complete separation 

is preferable. So, in contrast to the existing code, this type of experiments is now handled by

 a completely separated program package. Of course, many programs "overlap" or are even identical, but it nevertheless makes sense to treat the two kind of experiments (real-case versus idealised) separately.

\subsection{User Friendliness}

Computer programs are not written for a computer's joy, but because they should solve a problem

which cannot be solved analytically by any person. So, computer programs are tools which should 

make "life" easier for their users, or at least they should make some problems tractable. If the

user-interface of a computer program is chaotic and badly structured, the computer does not complain. 

The user, on the other hand, might give up using the tool simply because it needs too much effort to understand how the program must be handled. Scientific programs are especially in danger of 

having bad user-interfaces because these kind of programs are always "in the flow". This in turn

means that often it is not worthwhile to develop a graphical user-interface, as it is an absolute

must for commercial software. \\

\noindent

In the existing code for the PV inversion the user-interface was highly "chaotic". It consisted of 

several hundred lines of Linux Shell scripts which were difficult to read -partly due to lack of documentation, partly due to the unnecessarily complex structure. The behaviour of the inversion

was controlled by many variables set in the beginning of the Linux Shell script, set without

knowing exactly where and when these variables are used. This gave the user a "bad feeling" about

his experiments, simply because it was not always clear what changing parameters really meant.\\

\noindent

Although a graphical user-interface is not provided for the new version, an attractive command-based

interface is now available. This interface essentially consists of two parts. Firstly, a Linux Shell script

is provided which handles all calls to the computing Fortran programs. For instance, {\em inversion.sh prep} runs through all preparatory steps. User friendliness is reached by the fact the the user needs to make no changes to this Linux Shell script, although its limited size and documentation would make it quite simple.

Probably the best user-interface is provided by a parameter file. Such a file is provided in the new

version, and all parameters describing the inversion problem can be entered into this file. In fact, 

this kind of user-interface was motivated by sophisticated numerical weather prediction (NWP) models. There too,

the specifications characterising a model run are passed to the NWP model by means  of a parameter

file.\\

\noindent

What else made the existing user-interface chaotic? A close inspection revealed that the flow of data and information was absolutely non-transparent. Meteorological fields were moved from one file to 

another, only to be renamed there and subsequently being moved again to another file. For a part-time user it would have been impossible to keep track of the flow of information. In the new code the flow of information is very strictly controlled (see also the above subsection~3.1). In fact, there are many different files involved in the PV inversion:

input files, several intermediate files and finally the output files. It is not the number of files which makes a

process difficult to understand. If the flow of information is well controlled, it remains tractable. For the case of

the PV inversion, this information flow can be kept clear, and this was indeed an important aspect of the re-coding.  

% ----------------------------------------------------------------------------------------------

% PV Inversion for Real Cases 

% ----------------------------------------------------------------------------------------------

\newpage

\section{Inversion of an Extra-tropical PV Streamer}

\subsection{Scientific question}

The PV inversion starts with the selection of a distinct PV anomaly which should be

removed or added to the original PV distribution. It is then the aim to analyse

and compare the meteorological fields of horizontal wind and potential temperature

associated with the modified PV distribution with the corresponding fields associated with

the original PV distribution.\\

\begin{figure}[t]

\begin{center}

\begin{minipage}[t]{7.5cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{pv_20060114_18_315K.eps}

\end{minipage}

\hspace{0.5cm}

\begin{minipage}[t]{7.5cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{pv_20060115_18_315K.eps}

\end{minipage}

\\

\begin{minipage}[t]{7.5cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{pv_20060116_18_315K.eps}

\end{minipage}

\hspace{0.5cm}

\begin{minipage}[t]{7.5cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{pv_20060117_18_315K.eps}

\end{minipage}

\\

\begin{minipage}[t]{11cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{pv_colorbar.eps}

\end{minipage}

\\

\caption{\it Time evolution of Ertel's potential vorticity (PV) on the 315\,K isentropic surface. The four

plots are for 14, 15, 16 and 17 January 2006, 18\,UTC. The PV inversion will be exemplified in this user guide for 16 January 2006, 18\,UTC.}

\label{question}

\end{center}

\end{figure}

\noindent

The time evolution of potential vorticity (PV) on the 315\,K isentropic surface shows a distinct feature over the Eastern United States (Fig.\,\ref{question}). This so-called stratospheric PV streamer is indicative for

a deep intrusion of stratospheric air towards the equator. In the course of the four days shown, the

PV streamer moves towards the east, and in the later stages exhibits an cyclonic rolling-up. An eminent

question is how this stratospheric PV streamer influences the atmospheric flow in its neighbourhood. From theoretical

considerations, it can be expected that the impact on the horizontal velocity and on the static stability

is quite large (see introduction, Fig.\,2). By means of the so-called quasi-geostrophic omega equation (Holton, 1992), it is also expected that pronounced

vertical motions are directly associated with the passage of the PV streamer. It is the aim of this study -and of PV inversion in general- to surmount the qualitative argumentation, and to 

quantitatively assess the impact of such a structure on the atmospheric flow. \\

\noindent 

In the following chapters, all steps are discussed which are needed to quantify this impact. Here, we

exemplarily consider the date of 16 January 2006, 18\,UTC, and determine explicitly for this date 

the atmospheric response the upper-level (near tropopause level) PV structure.

\subsection{Necessary input files}

\begin{figure}[t]

\vspace{-2cm}

\begin{center}

\begin{minipage}[t]{8.0cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{q_20060116_18_500hPa.eps}

\end{minipage}

\hspace{-0.6cm}

\begin{minipage}[t]{8.0cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{t_20060116_18_500hPa.eps}

\end{minipage}

\\[-3cm]

\begin{minipage}[t]{8.0cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{uv_20060116_18_500hPa.eps}

\end{minipage}

\hspace{-0.6cm}

\begin{minipage}[t]{8.0cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{ps_20060116_18_500hPa.eps}

\end{minipage}

\\[-2cm]

\caption{\it Specific humidity (g/kg), temperature (in $^\circ$C), wind vectors

at 500\,hPa and surface pressure for 16 January 2006, 18\,UTC. These field are needed

as input for the PV inversion, and are available on the so-called P file.}

\label{inpp}

\end{center}

\end{figure}

Several preparatory steps are necessary before the PV inversion itself can be

performed. Here each of these steps is described in detail and -where possible- 

illustrated with suitable figures. Starting point of the analysis is the P and

Z file of the ECMWF (re-)analysis:

\begin{itemize}

\item {\bf P20060116\_18} with temperature T (in $^\circ$C), horizontal wind components U,V (in m/s) in

zonal and meridional direction, respectively,

      specific humidity Q (in kg/kg) and surface pressure PS (in hPa). 

\item {\bf Z20060116\_18} with geopotential height Z (in m) on a stack of pressure 

      levels (for instance on 850, 500 and 250\,hPa).

\end{itemize}

\begin{figure}[t]

\vspace{-2cm}

\begin{center}

\begin{minipage}[t]{8.0cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{z_20060116_18_250hPa.eps}

\end{minipage}

\hspace{-0.6cm}

\begin{minipage}[t]{8.0cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{z_20060116_18_500hPa.eps}

\end{minipage}

\\[-3.5cm]

\begin{minipage}[t]{8.0cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{z_20060116_18_850hPa.eps}

\end{minipage}

\\[-1.5cm]

\caption{\it Geopotential height (in m) on 250, 500 and 850\,hPa  for 16 January 2006, 18\,UTC. These field are needed

as input for the PV inversion, and are available on the Z file.}

\label{inpz}

\end{center}

\end{figure}

\noindent

Figure\,\ref{inpp} shows temperature T, specific humidity Q and wind vectors (U,V) at 500\,hPa and the surface pressure for 16 January 2006, 18\,UTC. A narrow band of high specific humidity is discernible on the downstream

side (to the east) of the PV streamer (upper-left panel). Most probably, this extended band is directly associated with

the southerly winds which prevail to the east of the PV streamer, and which are able to transport moist

air from the warm subtropics to the north. In the temperature field (upper-right panel), the PV streamer's impact is also discernible.  It is associated -at this level- by a local warm anomaly. On theoretical grounds, we would expect a local cold anomaly below the streamer and a local warm anomaly within and above the streamer (see introduction, Fig.\,2). The present signal indicates that the streamer reaches far down into the troposphere, and enforces a local warm anomaly at the low level of 500\,hPa. As already mentioned before, and explicitly shown by the wind vectors (lower-left panel), the PV streamer is also associated with a cyclonic flow. Besides the impact on the temperature field, this PV induced flow field is the most eminent impact of a PV streamer. Indeed, it is the wind field where the far-field effect of a PV streamer becomes most evident, because the streamer's impact on temperature and stratification is essentially confined to the

regions just below or above. 

Finally, the P file contains the surface pressure PS (lower-right panel), which of course to the largest extent simply reflects the surface topography.  So for instance,

the high topography of the Rocky Mountains or of the Greenland ice shield, is associated with low

surface pressure, 700\,hPa corresponding to about 3\,km height above sea level.\\

\noindent

In addition to the P file, a so-called Z file must be provided as input for the PV inversion. This file contains the geopotential height on a distinct set of pressure levels. Typically, these levels are

850, 500 and 250\,hPa. A plot of these levels is shown in Fig.\,\ref{inpz}. The need for this reference levels

will become evident in a subsequent section. In short: They are needed to integrate the hydrostatic equation and thereby to transform from pressure to geometrical height as vertical co-ordinate.

Note how the PV streamer over the west Atlantic is discernible in the geopotential as a deep trough. 

\subsection{An overview of the inversion}

\begin{figure}[t]

\begin{center}

\begin{minipage}[t]{6cm}

  \ForceWidth{1.0\textwidth}

  \BoxedEPSF{threesteps.eps}

\end{minipage}

\\

\caption{\it Data flow diagram for the three steps of the PV inversion. The preparatory steps are

essentially performed in the input directory, the inversion in the run directory and the output is

finally written to the output directory.}

\label{threesteps}

\end{center}

\end{figure}

The PV inversion can be split into three well-separated steps. The first step is associated with

preparations, in particular with the definition of the modified Ertel-PV field. In addition, some 

preparatory steps have to be performed in order to adapt the ECMWF grid to the cartesian grid which is 

assumed in the inversion algorithm. In the second step the atmospheric state is iteratively adjusted to

the modified PV field, as it was specified in the previous preparatory step. This second part comprises the "core" of the algorithm, i.e. the numerical inversion of the quasi-geostrophic potential 

vorticity equation by means of a successive over-relaxation (SOR) technique. Finally, in the third step

the modified atmospheric state is brought back from the cartesian grid of the inversion to the

original latitude/longitude grid of the ECMWF, including its hybrid vertical co-ordinate. After this

step a direct comparison between input and output fields is feasible due to the equivalent

grid structure of input and output fields. The basic three steps are shown in Fig.\,\ref{threesteps}, and will subsequently be described in greater detail in the following sections. Before doing so, a note has to be made on the user-interface for the PV inversion.\\

\noindent

The inversion is controlled by one master Linux Shell script (inversion.sh), which calls the performing Fortran programs

(for a complete listing of the Linux Shell script consider Appendix~9.2). 

A major aim of this work was to make the 

application of the inversion as easy as possible. Therefore, the typical user is not forced to study

the master script itself, but can provide all needed parameters by means of a so-called "parameter file"

(inversion.param).

The contents of this file is parsed by a Perl script, which extracts all needed information. In the

description of the subsequent steps, the relevant parts of this parameter file will be reproduced and

described. A complete listing of the parameter file can be found in Appendix~9.3.\\ 

\noindent

Nevertheless, it will be necessary for an advanced usage of the inversion package to go into

some greater details. For instance, if particular PV structures should be removed or added, the

implemented simple filtering approach (to be discussed in section~4.4.4) might not be sufficient.

In this case, the advanced user is encouraged to modify the Fortran code in order to fulfill his particular

tasks.\\ 

\noindent

In a first step some directories must be specified. This is done in the parameter file {\em inversion.param} in the section DATA. For the example of this study this section looks like:\\[-0.3cm]

\begin{center}

\begin{minipage}[t]{13cm}

\ForceWidth{1.0\textwidth}

\begin{verbatim}

BEGIN DATA

  DATE    = 20060116_18;                                

  INP_DIR = /lhome/sprenger/PV_Inversion_Tool/real/inp;   

  RUN_DIR = /lhome/sprenger/PV_Inversion_Tool/real/run; 

  OUT_DIR = /lhome/sprenger/PV_Inversion_Tool/real/out;

END DATA