3 
michaesp 
1 
\documentclass[12pt]{article}



2 



3 
\renewcommand{\topfraction}{0.93}



4 
\renewcommand{\bottomfraction}{0.93}



5 
\renewcommand{\textfraction}{.07}



6 



7 
\input BoxedEPS



8 
\SetRokickiEPSFSpecial



9 
\ForceOn



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



11 
\HideDisplacementBoxes



12 



13 
\textwidth16.0cm



14 
\textheight26.0cm



15 
\oddsidemargin0.5cm



16 
\topmargin1.5cm



17 
\headsep0cm



18 
\topskip0cm



19 



20 
\renewcommand{\baselinestretch}{1.2}



21 



22 
\begin{document}



23 



24 
\pagenumbering{Roman}



25 



26 
% 



27 
% Title page



28 
% 



29 



30 
\thispagestyle{empty}



31 



32 
\vspace{0.7cm}



33 



34 
\begin{center}



35 



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



37 
A user guide for realcase experiments}



38 



39 
\vspace{1cm}



40 
\begin{center}



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



42 
\ForceWidth{1.0\textwidth}



43 
\BoxedEPSF{Fig05_pv_320K.eps}



44 
\end{minipage}



45 
\end{center}



46 



47 



48 
\vspace{1cm}



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



50 
\vspace{0.3cm}



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



52 
\vspace{1.8cm}



53 
\today \\



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



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



56 
\end{center}



57 
\vspace{2cm}



58 



59 
\noindent



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



61 
Michael Sprenger\\



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



63 
ETH Zurich \\



64 
CH8092 Zurich, Switzerland \\



65 
email: michael.sprenger@env.ethz.ch



66 



67 



68 
\newpage



69 
\vspace{20cm}



70 
\noindent



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



72 



73 
% 



74 
% Abstract



75 
% 



76 



77 
\newpage



78 



79 
\begin{center}



80 
{\large \bf Abstract}



81 
\end{center}



82 



83 
\noindent



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



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



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



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



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



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



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



91 
as the invertibility and partitioning principle, respectively.\\



92 



93 
\noindent



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



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



96 
this is done with a socalled PV inversion, which in turn comprises several different steps, as for



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



98 
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 socalled



99 
quasigeostrophic balance condition. More specifically, the linear quasigeostrophic PV equation is



100 
solved. Since quasigeostrophic and Ertel's PV do not exactly coincide, an iterative technique



101 
is adopted to approach the nonlinear ErtelPV inversion by means of successively quasigeostrophic



102 
inversions.\\



103 



104 
\noindent



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



106 
Special focus was given to a clear and userfriendly program package, where all steps are



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



108 
inversion is controlled by one single Linux Shell script.



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



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



111 



112 
\newpage



113 



114 
\begin{center}



115 
{\large \bf Zusammenfassung}



116 
\end{center}



117 



118 
\noindent



119 
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,



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



121 
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



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



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



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



125 
Partitionierungsprinzip.\\



126 



127 
\noindent



128 
In dieser Arbeit wird ein Programmpaket vorgestellt, mit dessen Hilfe sich einzelne PVElemente



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



130 
lassen. Programmtechnisch besteht diese sogenannte PVInversion aus mehreren Schritten. So muss



131 
zum Beispiel eine Koordinatentransformation durchgef\"uhrt werden und eine PoissonGleichung



132 
numerisch gel\"ost werden. Damit die PVInversion sinnvoll ist, m\"ussen zwei Probleme gel\"ost



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



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



135 
stellt. In dem Programmpaket dieser Arbeit wird die sogenannte quasigeostrophische Balance



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



137 
linearen quasigeostrophischen PVGleichung. Das Problem der Nichtlinearit\"at wird dadurch gel\"ost, dass



138 
die Berechnung des genannten linearen Inversionsproblems mehrmals wiederholt wird. Iterativ



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



140 



141 
\noindent



142 
Ziel dieser Arbeit ist eine detailierte Darstellung aller Schritte, die bei einer PVInversion



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



144 
Strukturierung und Anwenderfreundlichkeit gelegt. Dadurch kann das Paket einerseits als BlackBox



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



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



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



148 
meisten Parameter, welche das gestellte Problem der PVInversion beschreiben, sind in einer



149 
separaten Parameterdatei eingetragen.



150 



151 



152 



153 



154 
% 



155 
% Inhalt



156 
% 



157 



158 
\newpage



159 
\section*{Contents}



160 



161 
\vspace{0.5cm}



162 
\renewcommand{\contentsname}{}



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



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



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



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



167 
\vspace{1.5cm}



168 
\tableofcontents



169 



170 



171 



172 
% 



173 
% Abkürzungen



174 
% 



175 



176 
\newpage



177 
\section*{Abbreviations}



178 



179 
\vspace{0.5cm}



180 
\begin{tabular}{rll}



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



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



183 
{\bf [3]} & {\bf ERA40} & Reanalysis of the ECMWF for the periods 19582001\\[2mm]



184 
{\bf [4]} & {\bf ERA15} & Reanalysis of the ECMWF for the periods 19791993\\[2mm]



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



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



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



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



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



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



191 
{\bf [11]} & {\bf netcdf } & Compact and randomaccess file format\\[2mm]



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



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



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



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



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



197 
{\bf [17]} & {\bf SOR} & Successive overrelaxation\\[2mm]



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



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



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



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



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



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



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



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



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



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



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



209 
\end{tabular}



210 



211 
\newpage



212 
\listoffigures



213 



214 
% 



215 
% Introduction and Motivation



216 
% 



217 



218 
\newpage



219 
\clearpage



220 
\pagenumbering{arabic}



221 



222 
\section{Introduction and Motivation}



223 



224 
\begin{figure}[t]



225 
\begin{center}



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



227 
\ForceWidth{1.0\textwidth}



228 
\BoxedEPSF{sat_wv.ps}



229 
\end{minipage}



230 
\\



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



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



233 
indicative for a lot of moisture in the upper troposphere. Note that highPV regions are coincident with



234 
dry regions. }



235 
\label{satwv}



236 
\end{center}



237 
\end{figure}



238 



239 
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 coordinate, or geopotential height, if pressure is used as the vertical coordinate. 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.\\



240 



241 
\noindent



242 
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 divergentfree and twodimensional flow its specification is sufficient to deduce the complete velocity field. We could call this the invertibility principle for vorticity in such a nondivergent and twodimensional 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).\\



243 



244 
\noindent



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



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



247 
generalisation of the barotropicvorticity equation for the real atmosphere is by far not trivial, but would be highly desirable.



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



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



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



251 
in twodimensional flow was first stated in a rough way by Kleinschmidt (1950). He was able to attribute



252 
some lowlevel flow features to an upperlevel PV anomaly, in his words to a "Zyklonenk\"orper". Figure\,1



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



254 
which form an elongated and narrow filament of high PV extending from Scandinavia south to the northwestern



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



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



257 



258 
\begin{figure}[t]



259 
\begin{center}



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



261 
\ForceWidth{1.0\textwidth}



262 
\BoxedEPSF{Fig11_PVInv_tropo1.ps}



263 
\end{minipage}



264 
\\



265 
\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 RossbyErtel potential vorticity has a strong discontinuity, by a factor of 6. The bold line is the tropopause, the horizontal



266 
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,



267 
1997).}



268 
\label{inversion}



269 
\end{center}



270 
\end{figure}



271 



272 
\noindent



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



274 
processes, whilst neglecting smallerscale features. Although these equations nowadays are considered no



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



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



277 
Along with this set of equations came a new version of potential vorticity. This quasigeostrophic



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



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



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



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



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



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



284 
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 socalled isentropes) are



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



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



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



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



289 
associated with such upperlevel distortions of the tropopause, and hence with a PV anomaly. Moreover, note



290 
that the upperlevel PV anomaly is not only associated with a deformation of the isentropes, but is



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



292 
is cyclonically (anticlockwise) oriented. The wind field is strongest at the tropopause level, but is even discernible at



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



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



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



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



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



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



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



300 



301 
\begin{figure}[t]



302 
\begin{center}



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



304 
\vspace{2cm}



305 
\ForceWidth{1.0\textwidth}



306 
\BoxedEPSF{streamer.eps}



307 
\end{minipage}



308 
\\[2cm]



309 
\caption{\it Stratospheric PV streamer and positions with stratospheretroposphere mass exchange. The section is



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



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



312 
\label{streamer}



313 
\end{center}



314 
\end{figure}



315 



316 



317 
\noindent



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



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



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



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



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



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



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



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



326 
features are quite common in the extratropics and are important for the mass exchange between the stratosphere



327 
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



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



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



330 
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



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



332 



333 



334 
\noindent



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



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



337 
key aspects for the restructuring of the program package.



338 
Then, in chapter 4 a detailed



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



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



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



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



343 



344 



345 



346 



347 
\newpage



348 



349 
\section{Mathematical Framework and Numerical Aspects}



350 



351 
In this study we use the quasigeostrophic approximation for the realcase



352 
inversion problem. In this limit the relative (quasigeostrophic) potential vorticity



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



354 



355 
\[



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



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



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



359 
\]



360 



361 
\vspace{0.3cm}



362 
\noindent



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



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



365 



366 
\[



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



368 
\quad \quad



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



370 
\quad \quad



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



372 
\]



373 



374 
\vspace{0.3cm}



375 
\noindent



376 
Here $N_0^2$, $\rho_0$ and $\theta_0$ denote the squared BruntV\"as\"ala frequency, the density and the



377 
potential temperature of a reference state depending only on the vertical coordinate $z$. $\psi$ is the



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



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



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



381 



382 



383 



384 
\noindent



385 
The above equations constitute a socalled von Neumann boundary problem for the elliptic



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



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



388 



389 
\[



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



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



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



393 
\sigma = 0



394 
\]



395 



396 
\vspace{0.3cm}



397 
\noindent



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



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



400 



401 
\[



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



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



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



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



406 
\]



407 



408 
\vspace{0.3cm}



409 
\noindent



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



411 
fields (quasigeostrophic PV $q$, streamfunction $\psi$,...) are defined on a threedimensional grid, whose



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



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



414 



415 
\[



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



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



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



419 
S = 0



420 
\]



421 



422 
\vspace{0.3cm}



423 
\noindent



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



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



426 



427 
\begin{eqnarray*}



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



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



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



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



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



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



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



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



436 
C_{i,j,k1/2} \cdot (\psi_{i,j,k}\psi_{i,j,k1})



437 
\end{eqnarray*}



438 



439 
\vspace{0.3cm}



440 
\noindent



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



442 



443 
\begin{eqnarray*}



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



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



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



447 
\end{eqnarray*}



448 



449 
\begin{figure}[t]



450 
\begin{center}



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



452 
\ForceWidth{1.0\textwidth}



453 
\BoxedEPSF{numericalgrid.eps}



454 
\end{minipage}



455 
\\



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



457 
in the xdirection (and correspondingly $\Delta y$ in the ydirection), the vertical grid spacing $\Delta z$. Note that in the vertical, a staggered grid at halflevels is also needed (taken from Fehlmann, 1997).}



458 
\label{grid}



459 
\end{center}



460 
\end{figure}



461 



462 
\vspace{0.3cm}



463 
\noindent



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



465 



466 
\[



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



468 
\]



469 



470 
\vspace{0.3cm}



471 
\noindent



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



473 
not only on the 3d grid where the quasigeostrophic PV and the streamfunction are defined, but to



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



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



476 
quasigeostrophic PV and the streamfunction is defined. The right side gives the indices for the



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



478 



479 
\noindent



480 
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):



481 



482 
\begin{eqnarray*}



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



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



485 
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)



486 
\end{eqnarray*}



487 



488 
\vspace{0.3cm}



489 
\noindent



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



491 



492 
\[



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



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



495 
\]



496 



497 
\vspace{0.3cm}



498 
\noindent



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



500 
all these definitions the fully discretised quasigeostrophic PV equation can be written as:



501 



502 
\begin{eqnarray*}



503 
S_{i,j,k} & = &



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



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



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



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



508 
\end{eqnarray*}



509 



510 
\vspace{0.3cm}



511 
\noindent



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



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



514 
von Neumann boundary conditions are:



515 



516 
\begin{eqnarray*}



517 
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]



518 
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]



519 
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]



520 
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)



521 
\end{eqnarray*}



522 



523 
\vspace{0.3cm}



524 
\noindent



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



526 
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:



527 



528 
\begin{eqnarray*}



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



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



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



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



533 
\end{eqnarray*}



534 



535 
\vspace{0.3cm}



536 
\noindent



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



538 
lower and upper lid of the domain.\\



539 



540 
\noindent



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



542 



543 
\[



544 
B\psi = b



545 
\]



546 



547 
\vspace{0.3cm}



548 
\noindent



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



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



551 
system has a solution if the following condition is fulfilled:



552 



553 
\[



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



555 
\]



556 



557 
\vspace{0.3cm}



558 
\noindent



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



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



561 
onedimensional because the nontrivial 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 onedimensional, the image of the operator $B$ is



562 
at most (m1) dimensional. Moreover, the operator $B$ is normal, and therefore the image is orthogonal to the



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



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



565 
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



566 
forcing terms.\\



567 



568 
\noindent



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



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



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



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



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



574 



575 
\noindent



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



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



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



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



580 



581 
\[



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



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



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



585 
\]



586 



587 
\vspace{0.3cm}



588 
\noindent



589 
converge toward the solution of the system $A\psi + \psi =b$. If we choose $A=B1$, the iterations converge toward a solution of the quasigeostrophic PV equation. \\



590 



591 
\noindent



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



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



594 
{\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.



595 



596 
% 



597 
% Program philosophy



598 
% 



599 



600 
\newpage



601 
\section{ReStructuring of the Code}



602 



603 
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 preexisting package. So, the question arises what added value this redevelopment and restructuring of the code is associated with. The basic



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



605 



606 
\begin{center}



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



608 
\begin{itemize}



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



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



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



612 
\end{itemize}



613 
\end{minipage}



614 
\end{center}



615 



616 
\vspace{0.3cm}



617 
\noindent



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



619 
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



620 
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



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



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



623 
a complete recoding is easier to be done than a restructuring of the existing version. The



624 
remedy against this code degeneration is a continuous restructuring. Hence, if changes in the code



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



626 
shortsightedly. The overall structure of the program should be kept in mind, and if possible



627 
longterm perspectives in code maintenance should be allowed, although such a longterm perspective



628 
momentarily leads to an additional effort. \\



629 



630 
\noindent



631 
Computer science, in particular software engineering, has defined the term "restructuring" (or



632 
"refactoring" for objectoriented programming languages) for the process how code should be



633 
maintained in order to avoid severe degeneration. Fowler defines the term refactoring in the



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



635 



636 
\begin{center}



637 
\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."}



638 
\end{center}



639 



640 
\vspace{0.6cm}



641 
\noindent



642 
Hence, the key aspect of restructuring is to continuously maintain the code, not only in



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



644 
changes.\\



645 



646 
\vspace{0.3cm}



647 
\noindent



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



649 
to the fact that cuttingedge science is, by "definition", investigating what is not known.



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



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



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



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



654 
so unorganised. Several works have been done with the program, several researchers included



655 
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



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



657 



658 
\noindent



659 
A first rough inspection immediately made



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



661 
Would it have been possible to only reorganise the code and write a new user guide? In a first



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



663 
recoding turned out to be more suitable. Hence, the original software package joined the



664 
sad fate of so many degenerated codes: a complete recoding. On the other hand, such a recoding



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



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



667 



668 
\noindent



669 
In the following three



670 
subsections, some key aspects of this recoding will be presented. They are by no means



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



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



673 



674 



675 



676 
\subsection{Transparency and Modularity}



677 



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



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



680 
for the PV inversion a classic three step splitting is appropriate: preprocessing, PV inversion, and postprocessing. Moreover, each of these three primary steps can further be split into several secondary substeps.



681 
In the existing code, the splitting of the problem into distinct subproblems was not clearly discernible. Indeed, the main Fortran program included not only the numerical inversion of the quasigeostrophic PV equation, but also some of the



682 
preparatory steps and some of the postprocessing steps. \\



683 



684 
\noindent



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



686 
separation is supported by the fact that three separated directories are used: There is one directory where preprocessing is done, on directory where inversion is done, and finally one directory where postprocessing is done. A flow of data between the three directories is allowed only at welldefined



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



688 
are moved from the run directory to the postprocessing directory.\\



689 



690 
\noindent



691 
Additional improvement resulted from a very clear separation of the subprocesses in the three mainprocesses (preprocessing, inversion, postprocessing). In fact, it was not tried to incorporate



692 
all preparatory steps into one large Fortran program, but instead each welldefined task (as for example transformation into a new coordinate system) constitutes a separate Fortran program. This



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



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



695 



696 
\noindent



697 
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: incode documentation. Every



698 
small subsection of a computer program should be documented. It should be possible to gain



699 
insight into an algorithm only by looking at the incode documentation. It is definitely not



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



701 
quite often even this most basic principle is not fulfilled. Hence, in the recoding of the



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



703 



704 
\noindent



705 
Finally a remark concerning the used programming language. A complete recoding of a software package might also be



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



707 
Fortran 77, and this language was also used for the recoding. Fortran, and in particular its older version



708 
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



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



710 
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



711 
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 objectoriented concepts offered by C++ and other modern languages. The PV inversion should be available



712 
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 recoding was done in Fortran is its high performance. In fact, PV inversion is very resource demanding, and compilers must create fastrunning codes in order



713 
the PV inversion tool to be helpful.\\



714 



715 
\noindent



716 
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 wellstructured and carefully developed Fortran program is



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



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



719 
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



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



721 
code lacks clarity due to the large number of subroutines. In the present recoding of the PV inversion, focus was given



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



723 



724 



725 
\subsection{Universality and Model Independence}



726 



727 
Universality was another aspect which was important in recoding the PV inversion. In the present



728 
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,



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



730 
model (EM) of the German weather service (DWD). Later it was applied to the higherresolution



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



732 
nonhydrostatic local model (LM) developed by the DWD replaced the HRM in regional



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



734 
Moreover, PV inversion is particularly attractive for globalscale data sets. So, the method and the



735 
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



736 
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 terrainfollowing coordinate system which in higher levels becomes a pressurelevel coordinate system, quite in contrast what is used by the regional LM model (a variant of geometrical height). \\



737 



738 
\noindent



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



740 
dealt with. In the recoding, model aspects were removed in the first two preparatory steps (out of



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



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



743 
transforms the fields into a local quasicartesian coordinate 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



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



745 
problem by introducing a new quasicartesian coordinate system centered over the North Pole. In this



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



747 



748 
\noindent



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



750 
is generally so different from what is needed for realcase studies that a complete separation



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



752 
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 (realcase versus idealised) separately.



753 



754 



755 



756 



757 



758 
\subsection{User Friendliness}



759 



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



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



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



763 
userinterface of a computer program is chaotic and badly structured, the computer does not complain.



764 
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



765 
having bad userinterfaces because these kind of programs are always "in the flow". This in turn



766 
means that often it is not worthwhile to develop a graphical userinterface, as it is an absolute



767 
must for commercial software. \\



768 



769 
\noindent



770 
In the existing code for the PV inversion the userinterface was highly "chaotic". It consisted of



771 
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



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



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



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



775 



776 
\noindent



777 
Although a graphical userinterface is not provided for the new version, an attractive commandbased



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



779 
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.



780 
Probably the best userinterface is provided by a parameter file. Such a file is provided in the new



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



782 
this kind of userinterface was motivated by sophisticated numerical weather prediction (NWP) models. There too,



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



784 
file.\\



785 



786 



787 
\noindent



788 
What else made the existing userinterface chaotic? A close inspection revealed that the flow of data and information was absolutely nontransparent. Meteorological fields were moved from one file to



789 
another, only to be renamed there and subsequently being moved again to another file. For a parttime 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:



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



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



792 
the PV inversion, this information flow can be kept clear, and this was indeed an important aspect of the recoding.



793 



794 



795 



796 



797 
% 



798 
% PV Inversion for Real Cases



799 
% 



800 



801 
\newpage



802 



803 
\section{Inversion of an Extratropical PV Streamer}



804 



805 
\subsection{Scientific question}



806 



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



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



809 
and compare the meteorological fields of horizontal wind and potential temperature



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



811 
the original PV distribution.\\



812 



813 
\begin{figure}[t]



814 
\begin{center}



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



816 
\ForceWidth{1.0\textwidth}



817 
\BoxedEPSF{pv_20060114_18_315K.eps}



818 
\end{minipage}



819 
\hspace{0.5cm}



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



821 
\ForceWidth{1.0\textwidth}



822 
\BoxedEPSF{pv_20060115_18_315K.eps}



823 
\end{minipage}



824 
\\



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



826 
\ForceWidth{1.0\textwidth}



827 
\BoxedEPSF{pv_20060116_18_315K.eps}



828 
\end{minipage}



829 
\hspace{0.5cm}



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



831 
\ForceWidth{1.0\textwidth}



832 
\BoxedEPSF{pv_20060117_18_315K.eps}



833 
\end{minipage}



834 
\\



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



836 
\ForceWidth{1.0\textwidth}



837 
\BoxedEPSF{pv_colorbar.eps}



838 
\end{minipage}



839 
\\



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



841 
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.}



842 
\label{question}



843 
\end{center}



844 
\end{figure}



845 



846 
\noindent



847 
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 socalled stratospheric PV streamer is indicative for



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



849 
PV streamer moves towards the east, and in the later stages exhibits an cyclonic rollingup. An eminent



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



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



852 
is quite large (see introduction, Fig.\,2). By means of the socalled quasigeostrophic omega equation (Holton, 1992), it is also expected that pronounced



853 
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



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



855 



856 
\noindent



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



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



859 
the atmospheric response the upperlevel (near tropopause level) PV structure.



860 



861 



862 



863 
\subsection{Necessary input files}



864 



865 
\begin{figure}[t]



866 
\vspace{2cm}



867 
\begin{center}



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



869 
\ForceWidth{1.0\textwidth}



870 
\BoxedEPSF{q_20060116_18_500hPa.eps}



871 
\end{minipage}



872 
\hspace{0.6cm}



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



874 
\ForceWidth{1.0\textwidth}



875 
\BoxedEPSF{t_20060116_18_500hPa.eps}



876 
\end{minipage}



877 
\\[3cm]



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



879 
\ForceWidth{1.0\textwidth}



880 
\BoxedEPSF{uv_20060116_18_500hPa.eps}



881 
\end{minipage}



882 
\hspace{0.6cm}



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



884 
\ForceWidth{1.0\textwidth}



885 
\BoxedEPSF{ps_20060116_18_500hPa.eps}



886 
\end{minipage}



887 
\\[2cm]



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



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



890 
as input for the PV inversion, and are available on the socalled P file.}



891 
\label{inpp}



892 
\end{center}



893 
\end{figure}



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



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



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



897 
Z file of the ECMWF (re)analysis:



898 



899 
\begin{itemize}



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



901 
zonal and meridional direction, respectively,



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



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



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



905 
\end{itemize}



906 



907 
\begin{figure}[t]



908 
\vspace{2cm}



909 
\begin{center}



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



911 
\ForceWidth{1.0\textwidth}



912 
\BoxedEPSF{z_20060116_18_250hPa.eps}



913 
\end{minipage}



914 
\hspace{0.6cm}



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



916 
\ForceWidth{1.0\textwidth}



917 
\BoxedEPSF{z_20060116_18_500hPa.eps}



918 
\end{minipage}



919 
\\[3.5cm]



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



921 
\ForceWidth{1.0\textwidth}



922 
\BoxedEPSF{z_20060116_18_850hPa.eps}



923 
\end{minipage}



924 
\\[1.5cm]



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



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



927 
\label{inpz}



928 
\end{center}



929 
\end{figure}



930 



931 
\noindent



932 
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



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



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



935 
air from the warm subtropics to the north. In the temperature field (upperright 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 (lowerleft 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 farfield effect of a PV streamer becomes most evident, because the streamer's impact on temperature and stratification is essentially confined to the



936 
regions just below or above.



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



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



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



940 



941 
\noindent



942 
In addition to the P file, a socalled 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



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



944 
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 coordinate.



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



946 



947 



948 



949 
\subsection{An overview of the inversion}



950 



951 
\begin{figure}[t]



952 
\begin{center}



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



954 
\ForceWidth{1.0\textwidth}



955 
\BoxedEPSF{threesteps.eps}



956 
\end{minipage}



957 
\\



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



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



960 
finally written to the output directory.}



961 
\label{threesteps}



962 
\end{center}



963 
\end{figure}



964 



965 
The PV inversion can be split into three wellseparated steps. The first step is associated with



966 
preparations, in particular with the definition of the modified ErtelPV field. In addition, some



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



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



969 
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 quasigeostrophic potential



970 
vorticity equation by means of a successive overrelaxation (SOR) technique. Finally, in the third step



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



972 
original latitude/longitude grid of the ECMWF, including its hybrid vertical coordinate. After this



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



974 
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 userinterface for the PV inversion.\\



975 



976 



977 
\noindent



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



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



980 
A major aim of this work was to make the



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



982 
the master script itself, but can provide all needed parameters by means of a socalled "parameter file"



983 
(inversion.param).



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



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



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



987 



988 



989 
\noindent



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



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



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



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



994 
tasks.\\



995 



996 



997 
\noindent



998 
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]



999 
\begin{center}



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



1001 
\ForceWidth{1.0\textwidth}



1002 
\begin{verbatim}



1003 
BEGIN DATA



1004 
DATE = 20060116_18;



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



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



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



1008 
END DATA



1009 
\end{verbatim}



1010 
\end{minipage}



1011 
\end{center}



1012 



1013 
\vspace{0.3cm}



1014 
\noindent



1015 
The first argument (DATE) gives the date of the case. Then the following three entries specify the input (INP\_DIR), the run (RUN\_DIR) and the output (OUT\_DIR) directory, respectively. These directories can be created by the call



1016 
\\[0.3cm]



1017 
\begin{center}



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



1019 
\ForceWidth{1.0\textwidth}



1020 
\begin{verbatim}



1021 
inversion.sh inst



1022 
\end{verbatim}



1023 
\end{minipage}



1024 
\end{center}



1025 



1026 
\vspace{0.3cm}



1027 
\noindent



1028 
For test reasons a set of sample files is provided. These sample P20060115\_18 and Z20060115\_18 files can be copied to the input directory with the call



1029 
\\[0.3cm]



1030 
\begin{center}



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



1032 
\ForceWidth{1.0\textwidth}



1033 
\begin{verbatim}



1034 
inversion.sh sample



1035 
\end{verbatim}



1036 
\end{minipage}



1037 
\end{center}



1038 



1039 
\vspace{0.3cm}



1040 
\noindent



1041 
The fields available on the two files are:



1042 



1043 
\begin{center}



1044 
\begin{tabular}{lll}



1045 
\hline



1046 
Q & Specific humidity ($kg/kg$) & P20060116\_18 \\



1047 
U & Zonal velocity ($m/s$) & P20060116\_18 \\



1048 
V & Meridional velocity ($m/s$) & P20060116\_18 \\



1049 
T & Temperature ($^\circ C$) & P20060116\_18 \\



1050 
OMEGA & Vertical velocity ($Pa/s$) & P20060116\_18 \\



1051 
PS & Surface pressure ($hPa$) & P20060116\_18 \\



1052 
\hline



1053 
Z & Geopotential height ($m$) on 850, 500 and 250\,hPa & Z20060115\_18 \\



1054 
\hline



1055 
\end{tabular}



1056 
\end{center}



1057 



1058 



1059 
\vspace{0.5cm}



1060 
\noindent



1061 
With these preparatory steps, the inversion can be performed. The following sections discuss the three basic step, i.e. the preparation of the input files for the inversion, the inversion itself and the post processing of the output files.



1062 



1063 



1064 



1065 
\subsection{Preprocessing: Defining the inversion problem [prep]}



1066 



1067 
\begin{figure}[t]



1068 
\begin{center}



1069 
\begin{minipage}[t]{12cm}



1070 
\vspace{2.5cm}



1071 
\ForceWidth{1.0\textwidth}



1072 
\BoxedEPSF{preparation.eps}



1073 
\end{minipage}



1074 
\\



1075 
\caption{\it Flowchart of preparatory steps. Additionally, to the left of the single steps the



1076 
essential flow of information is shown. For instance, H $\rightarrow$ R means that the input is taken from



1077 
the file with prefix H and the output is then written to the file with prefix R.}



1078 
\label{preparation}



1079 
\end{center}



1080 
\end{figure}



1081 



1082 
Firstly, the vertical coordinate



1083 
has to be changed from pressure to geometrical height (section~4.4.1) and the fields have to be



1084 
transformed from the geographical latitude/longitude grid to a cartesian grid (section~4.4.2). Then additional meteorological



1085 
fields have to be computed on this new cartesian grid (section~4.4.3). With these two preparations the



1086 
ErtelPV anomaly can be defined (section~4.4.4), and the input files for the quasigeostrophic PV inversion be written



1087 
(section~4.4.5). Some additional steps are then: definition of a reference profile (in section~4.4.6) and



1088 
the transformation of the coastlines into the cartesian grid (section~4.4.6).\\



1089 



1090 
\noindent



1091 
These single steps have to be done once for every inversion problem. A complete flowchart of



1092 
the preparatory steps is shown in Fig.\,\ref{preparation}. In the following description every single step will



1093 
be launched and discussed individually. If desired, the call



1094 
\\[0.0cm]



1095 
\begin{center}



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



1097 
\ForceWidth{1.0\textwidth}



1098 
\begin{verbatim}



1099 
inversion.sh prep



1100 
\end{verbatim}



1101 
\end{minipage}



1102 
\end{center}



1103 



1104 
\vspace{0.6cm}



1105 
\noindent



1106 
will run through all steps without waiting for intermediate check. This is particularly practical if all



1107 
settings remain essentially unchanged, and therefore a detailed check after each step is not necessary.



1108 



1109 



1110 



1111 



1112 
\subsubsection{Transformation to height levels [prep1]}



1113 



1114 
\begin{figure}[t]



1115 
\begin{center}



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



1117 
\ForceWidth{1.0\textwidth}



1118 
\BoxedEPSF{sigma_level.eps}



1119 
\end{minipage}



1120 
\\



1121 
\caption{\it Illustration of the terrainfollowing hybrid coordinate system of the ECMWF grid.



1122 
Here it is shown with 15 levels, recent versions of the ECMWF model contain 60 hybrid levels for the ERA40 reanalysis and 91 levels for the operational analysis and deterministic 10\,day forecast [illustration taken



1123 
from "An Introduction to dynamic meteorology" by Holton, 1992].}



1124 
\label{sigma}



1125 
\end{center}



1126 
\end{figure}



1127 



1128 



1129 
The input data are available on a latitude/longitude grid



1130 
in the horizontal and on a socalled hybrid $\sigma$grid on the vertical. Whereas the former is



1131 
well known, the latter needs some explanation. The basic idea is illustrated in Fig.\,\ref{sigma}. The



1132 
model levels are terrainfollowing, and the pressure p(i,j,k) at a specific grid point (with grid indices i,j,k) is given by the expression:



1133 
\[



1134 
p(i,j,k) = ak(k) + bk(k)*ps(i,j)



1135 
\]



1136 



1137 
\vspace{0.3cm}



1138 
\noindent



1139 
where $ps(i,j)$ is the pressure at the surface at the specified grid position. The two onedimensional



1140 
arrays $ak$ and $bk$ are part of the P file, more specifically of the constants file associated with the P file. The main advantage of this hybrid grid is its nonintersection with the ground, which



1141 
facilitates many numerical calculations.\\



1142 



1143 
\noindent



1144 
The ECMWF fields are easily plotted on pressure surfaces since the vertical coordinate of the



1145 
underlying data set is essentially based upon pressure. On the other hand, the inversion program



1146 
assumes a cartesian grid with geometrical height (in m) as vertical coordinate. Therefore, all needed fields from the



1147 
original P file must be interpolated onto a stack of height levels. This transformation is easily



1148 
done by means of the hydrostatic equation:



1149 
\[



1150 
\frac{\partial p}{\partial z} =  \rho \cdot g



1151 
\]



1152 



1153 
\vspace{0.3cm}



1154 
\noindent



1155 
where $\rho$ denotes the density of moist air and $g$ is the Earth's gravity. Note that this



1156 
equation can be reformulated to one including temperature if the ideal gas equation



1157 
\[



1158 
p = \rho \cdot R_d \cdot T_v



1159 
\]



1160 
is used and $\rho$ replaced. Here $R_d$ is the ideal gas constant of dry air and $T_v$ is the virtual temperature,



1161 
i.e. the temperature corrected for the water vapour contents of the air. The definition of



1162 
virtual temperature takes into account the specific humidity of an air parcel and its (sensible)



1163 
temperature (see for instance Wallace and Hobbs, 1977):\\



1164 
\[



1165 
T_v = T \cdot \{ 1  \frac{q}{\epsilon} \cdot (1  \epsilon) \}^{1} \approx



1166 
T \cdot ( 1 + \epsilon \cdot q )



1167 
\]



1168 



1169 
\vspace{0.3cm}



1170 
\noindent



1171 
where $T$ is the temperature (in K), $q$ the specific humidity (in kg/kg), $\epsilon=R_d/R_v$ is the ratio of the gas constant for dry ($R_d$) air and water vapour ($R_v$). \\



1172 



1173 
\noindent



1174 
With the previous definitions, the hydrostatic equation can be integrated vertically to get



1175 
the height as a function of the pressure p:



1176 



1177 
\[



1178 
z(p) = z_{ref}  \frac{R_d}{g} \cdot \int_{\log(p_{ref})}^{\log(p)} T_v \cdot d\log(p)



1179 
\]



1180 



1181 
\vspace{0.3cm}



1182 
\noindent



1183 
where $z_{ref}$ and $p_{ref}$ denote the reference level. If several reference levels are given



1184 
(for example at 850, 500 and 250\,hPa, see Fig.\,7) the integration can be performed from the nearest



1185 
reference level to the specified pressure level. Note that at least one reference level



1186 
must be given, i.e. at one pressure level the geometrical height must be known.



1187 
Typically, this correspondence is known for the 500\,hPa level (see Fig.\,\ref{hydro}). As a result



1188 
of the integration, the geometrical height is given as a function of the pressure at all points



1189 
of the model grid.\\



1190 



1191 
\begin{figure}[t]



1192 
\begin{center}



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



1194 
\ForceWidth{1.0\textwidth}



1195 
\BoxedEPSF{pro_t_70w40n_20060115_18.eps}



1196 
\end{minipage}



1197 
\hspace{0.5cm}



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



1199 
\ForceWidth{1.0\textwidth}



1200 
\BoxedEPSF{pro_q_70w40n_20060115_18.eps}



1201 
\end{minipage}



1202 
\\



1203 
\caption{\it Left: Vertical profile of temperature (in $^\circ$C) in dependence of pressure at 70$^\circ$W



1204 
and 40$^\circ$ for 15 January 2006, 18\,UTC. Additionally the reference pressure levels are included with



1205 
the corresponding geopotential height (in m). Right: Vertical profile of specific humidity (in g/kg).}



1206 
\label{hydro}



1207 
\end{center}



1208 
\end{figure}



1209 



1210 
\noindent



1211 
In addition to the geometrical height at each grid point, the hydrostatic equation can be integrated



1212 
downward to obtain a consistent estimate of topography height. In fact, if the surface pressure



1213 
$p_s$ is given, the height of the topography is readily determined by\\



1214 



1215 
\[



1216 
z(p_s) = z_{ref}  \frac{R_d}{g} \cdot \int_{\log(p_{ref})}^{\log(p_s)} T_v \cdot d\log(p)



1217 
\]



1218 



1219 
\vspace{0.3cm}



1220 
\noindent



1221 
The advantage of the thus determined topography (instead of using the ECMWF topography) is its



1222 
consistency.\\



1223 



1224 
\noindent



1225 
Having determined the geometrical height $z(i,j,k)$ at each grid point $i,j,k$ in addition to the already known



1226 
pressure $p(i,j,k)$, it is straightforward to interpolate a field (temperature, velocity,...) onto an arbitrary height level.



1227 
The method adopted here is to lie a natural cubic spline through the grid points along a vertical



1228 
profile and then to evaluate the spline at the prespecified height levels. In detail: The cubic spline at a horizontal grid position $i,j$ is based upon the values $[x_k=z(i,j,k),y_k=T(i,j,k)]$ if temperature $T$ is to be interpolated onto



1229 
a stack of height levels. Having determined the associated cubic spline $T_{sp}$, temperature at an arbitrary height $z$



1230 
can be calculated by evaluating the cubic spline at this level: $T_{sp}(z)$. The algorithm for the cubic spline is taken from Press et al. (1992).



1231 
\\



1232 



1233 



1234 
\noindent



1235 
The interpolation onto height levels is done with the call\\[0.3cm]



1236 
\begin{center}



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



1238 
\begin{verbatim}



1239 
inversion.sh prep1



1240 
\end{verbatim}



1241 
\end{minipage}



1242 
\end{center}



1243 



1244 
\vspace{0.3cm}



1245 
\noindent



1246 
which numerically



1247 
integrates the hydrostatic equation and therefore attributes to each grid point not only its



1248 
pressure (in hPa) but also its geometrical height (in m). The relevant numerical parameters are taken from the



1249 
parameter file {\em inversion.param}, where the following section is essential:\\



1250 
\begin{center}



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



1252 
\begin{verbatim}



1253 
BEGIN GRID



1254 
GEO_ZMIN = 0.;



1255 
GEO_NZ = 125 ;



1256 
GEO_DZ = 200.;



1257 
END GRID



1258 
\end{verbatim}



1259 
\end{minipage}



1260 
\end{center}



1261 



1262 
\vspace{0.3cm}



1263 
\noindent



1264 
The three parameters describe the new vertical coordinate. The lowest levels is assumed to be at ground



1265 
(GEO\_ZMIN), the number of vertical levels is 125 (GEO\_NZ) and the vertical resolution is equidistantly



1266 
200\,m (GEO\_DZ). Hence, the new vertical grid spans the range from ground up to 25\,km.\\



1267 



1268 
\noindent



1269 
The output file H20060115\_18 is a new netcdf file which has now geometrical height as vertical coordinate



1270 
and the variables ({U,V,T,Q}) interpolated onto the height levels from surface up to 25\,km. Note that the height of the topography ({ORO}) is also



1271 
determined. If a model level is below the topography, all field values are set to missing data. The following table gives all the variables which are available on the new H file:\\



1272 



1273 
\begin{center}



1274 
\begin{tabular}{lll}



1275 
\hline



1276 
Q & Specific humidity ($kg/kg$) & H20060116\_18 \\



1277 
U & Zonal velocity ($m/s$) & H20060116\_18 \\



1278 
V & Meridional velocity ($m/s$) & H20060116\_18 \\



1279 
T & Temperature ($^\circ C$) & H20060116\_18 \\



1280 
P & Pressure ($hPa$) & H20060116\_18 \\



1281 
ORO & Height of topography ($m$) & H20060116\_18 \\



1282 
\hline



1283 
\end{tabular}



1284 
\end{center}



1285 



1286 
\begin{figure}[t]



1287 
\begin{center}



1288 
\begin{minipage}[t]{14cm}



1289 
\vspace{1cm}



1290 
\ForceWidth{1.0\textwidth}



1291 
\BoxedEPSF{latlon_nh.eps}



1292 
\end{minipage}



1293 
\\[0.5cm]



1294 
\caption{\it Latitude/longitude grid of the ECMWF P and Z files. For clarity, the grid spacing was reduced to 5 degrees, instead of the 1 degree in the files. Meteorological fields as for instance



1295 
temperature, wind vectors, pressure and geopotential are given at the intersection of the latitude and



1296 
longitude circles. Note how this grid becomes singular toward the north pole.}



1297 
\label{latlon}



1298 
\end{center}



1299 
\end{figure}



1300 



1301 



1302 
\subsubsection{Rotating to a quasicartesian coordinate system [prep2]}



1303 



1304 
\begin{figure}[t]



1305 
\begin{center}



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



1307 
\ForceWidth{1.0\textwidth}



1308 
\BoxedEPSF{local_cart_grid.eps}



1309 
\end{minipage}



1310 
\\



1311 
\caption{\it Local cartesian grid which is centered at the PV streamer's position and which



1312 
shows much less distortion due to the sphericity of the Earth. The center of the new grid



1313 
corresponds to the rotated latitude and longitude (0,0). The boundaries are at $31^\circ$ and



1314 
$+31^\circ$ rotated longitude, and at $31^\circ$ and $+31^\circ$ rotated latitude.}



1315 
\label{localcart}



1316 
\end{center}



1317 
\end{figure}



1318 



1319 
The input fields are given on a latitude/longitude grid, with after the previous step the fields



1320 
on an equidistant grid of height levels. The inversion itself assumes a cartesian grid, hence a



1321 
uniform grid spacing is also assumed in the horizontal directions. This is clearly not the case for the



1322 
latitude/longitude grid, as illustrated in Fig.\,\ref{latlon}.



1323 
Due to the sphericity of the Earth, the zonal



1324 
grid spacing decreases with increasing geographical latitude. As a remedy to this singularity, all



1325 
input fields are transformed to a "quasicartesian" grid with its center on the PV structure of



1326 
interest. The transformation between the two grids is illustrated in Fig.\,\ref{localcart}. It becomes evident



1327 
that the grid distortion is significantly reduced by this step. The formula which relates the two



1328 
grids is based upon the coordinate transformation given in the Appendix~A (Fortran subroutines {\em lmstolm} and



1329 
{\em phstoph}). The geographical latitude/longitude corresponding to a grid point in the rotated (quasicartesian) coordinate system is obtained in the following way. Firstly, the rotated latitude/longitude $\phi_{rot},\lambda_{rot}$ is transformed into



1330 



1331 
\begin{eqnarray*}



1332 
\lambda'_{rot} & = & 90+lmstolm(\phi_{rot},\lambda_{rot}90,90+\alpha,180) \\[0.2cm]



1333 
\phi'_{rot} & = & phstoph(\phi_{rot},\lambda_{rot}90,90+\alpha,180)



1334 
\end{eqnarray*}



1335 



1336 
\vspace{0.3cm}



1337 
\noindent



1338 
where $\alpha$ is the rotation angle given in the parameter file (CROT, see below). These two



1339 
values are then, in a second rotation, transformed into geographical latitude and longitude:



1340 



1341 
\begin{eqnarray*}



1342 
\lambda_{geo} & = & lmstolm(\phi'_{rot},\lambda'_{rot},90\phi_{cen},\lambda_{cen}180) \\[0.2cm]



1343 
\phi_{geo} & = & phstoph(\phi'_{rot},\lambda'_{rot},90\phi_{cen},\lambda_{cen}180))



1344 
\end{eqnarray*}



1345 



1346 
\vspace{0.3cm}



1347 
\noindent



1348 
where the centre of the PV anomaly (the centre of the quasicartesian coordinate system in geogrpahical



1349 
latitude/longitude coordinates) is given by $\phi_{cen}$ and $\lambda_{cen}$ (CLAT and CLON in the



1350 
parameter file, see below).\\



1351 



1352 
\begin{figure}[t]



1353 
\begin{center}



1354 
\begin{minipage}[t]{8cm}



1355 
\ForceWidth{1.0\textwidth}



1356 
\BoxedEPSF{latlon_grid.eps}



1357 
\end{minipage}



1358 
\hspace{0.5cm}



1359 
\begin{minipage}[t]{8cm}



1360 
\ForceWidth{1.0\textwidth}



1361 
\BoxedEPSF{xy_grid.eps}



1362 
\end{minipage}



1363 
\\



1364 
\caption{\it Quasicartesian grid which is centered at the PV streamer's position and which



1365 
shows much less distortion due to the sphericity of the Earth. On the left panel the latitude and



1366 
longitude on the Earth surface is shown, on the right panel the x and y coordinates in the new



1367 
coordinate system. Note how the distortion is considerably reduced with the new x,y coordinates (right), as compared with the original latitude/longitude coordinates (left). Additionally, the topography (in m above sea level) is shown



1368 
in color: Clearly discernible is Greenland and the east coast of North America.}



1369 
\label{xylatlon}



1370 
\end{center}



1371 
\end{figure}



1372 



1373 
\noindent



1374 
The transformation looks quite complicated, but has an easy geometrical interpretation. In fact, we create a new



1375 
latitude/longitude grid, but now longitude and latitude are defined not relative to the Earth's north



1376 
pole, but relative to a rotated imaginary "north pole". Therefore, we call the new latitudes/longitudes rotated coordinates. The link between the geographical latitude/longitude and



1377 
the rotated ones is shown in Fig.\,\ref{xylatlon} for the case study. Here, the left panel gives the latitude and longitude at the grid positions of the quasicartesian grid. The right panel gives the quasicartesian coordinates x and y, which correspond to the grandcircle distance on the Earth surface.



1378 
While the transformation of a scalar quantity (like for instance temperature, pressure, wind speed) is straightforward with the above formulas, the transformation of a vector field is much more



1379 
complicated. This is easily illustrated with the following example: Consider a perfectly zonal flow



1380 
in the geographical grid, i.e. there is only a wind component along the geographical latitude circles, but none along the geographical longitude circles (hence, in the left panel of Fig.\,14 the velocity vector is parallel to the latitude



1381 
circles). In the rotated grid the quasicartesian grid this



1382 
pure zonal velocity is transformed into a wind vector which has components both along the rotated



1383 
longitude and latitude circles (i.e. in Fig.\,14 the vector has a component along the horizontal and the



1384 
vertical axis).



1385 
In principle, a complicated formula could be derived which gives an explicit expression how the



1386 
individual velocity components transform. Here, we adopt a more pragmatic approach: We numerically



1387 
determine the local rotation angle between the two orthogonal coordinate systems, and subsequently



1388 
use these angles to get the transformed velocity components. \\



1389 



1390 
\begin{figure}[t]



1391 
\begin{center}



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



1393 
\ForceWidth{1.0\textwidth}



1394 
\BoxedEPSF{coriol_grid.eps}



1395 
\end{minipage}



1396 
\\



1397 
\caption{\it Coriolis parameter (in $10^{4}s^{1}$) in the quasicartesian coordinate frame.



1398 
The Coriolis parameter is a local measure for the Earth's rotation. It is constant along



1399 
latitude circles. Note that $f$ becomes maximal at the pole and vanishes at the equator. No



1400 
PV inversion can be performed if the cartesian grid crosses the equator and hence crosses the



1401 
zero isoline of the Coriolis parameter.}



1402 
\label{coriol}



1403 
\end{center}



1404 
\end{figure}



1405 



1406 
\noindent



1407 
Of course, the Coriolis parameter must be kept in the transformation to the quasicartesian



1408 
coordinate frame. Tis parameter is given by:



1409 



1410 
\[



1411 
f = 2 \cdot \Omega \cdot sin(\phi)



1412 
\]



1413 



1414 
\vspace{0.3cm}



1415 
\noindent



1416 
where $\Omega$ is the angular velocity of the Earth rotation and $\phi$ is the geographical latitude.



1417 
From the definition it is clear that the Coriolis parameter is constant along latitude circles. Furthermore,



1418 
it takes its maximum value at the pole and vanishes at the equator. A typical value for the midlatitudes



1419 
is $f=10^{4}s^{1}$. For the case study the Coriolis parameter is shown in Fig.\,\ref{coriol}.\\



1420 



1421 



1422 



1423 
\noindent



1424 
Numerically, the rotation of the fields to the quasicartesian coordinate frame is done with the



1425 
call\\[0.3cm]



1426 
\begin{center}



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



1428 
\begin{verbatim}



1429 
inversion.sh prep2



1430 
\end{verbatim}



1431 
\end{minipage}



1432 
\end{center}



1433 



1434 
\vspace{0.3cm}



1435 
\noindent



1436 
The program needs some information from the parameter file. The relevant section is:\\[0.3cm]



1437 
\begin{center}



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



1439 
\begin{verbatim}



1440 
BEGIN GRID



1441 
ROT_NX = 250 ;



1442 
ROT_NY = 250 ;



1443 
ROT_DX = 0.25;



1444 
ROT_DY = 0.25;



1445 
CLON = 65.;



1446 
CLAT = 45.;



1447 
CROT = 0.;



1448 
END GRID



1449 
\end{verbatim}



1450 
\end{minipage}



1451 
\end{center}



1452 



1453 
\vspace{0.3cm}



1454 
\noindent



1455 
The first two parameters (ROT\_NX and ROT\_NY) give the number of grid points in the horizontal directions for



1456 
the rotated coordinate system. Here it is assumed that the grid has 250 grid points in x and in y direction. Note



1457 
that the number of grid points in the z direction needs not to be specified because this information is



1458 
already available on the file H20060115\_18. The new horizontal resolution of the rotated grid is given by the



1459 
next two values. They give the resolution in x (ROT\_DX) and in y (ROT\_DY) direction. In the present case these



1460 
resolutions are 0.25 degrees, which corresponds to approximately 28\,km, according to the formula\\



1461 



1462 
\[



1463 
0.25 \cdot \frac{2\pi}{360} \cdot R_E \quad \quad \mbox{with the Earth's radius} \quad R_E=6370\,km



1464 
\]



1465 



1466 
\vspace{0.3cm}



1467 
\noindent



1468 
Note that the resolution of the rotated grid is higher than the resolution of the input grid, where a value of



1469 
1\,degree is given. Finally, the last three parameters specify where the new coordinate system is centered on the



1470 
globe. Here, this position is at 65\,W (CLON) and 45\,N (CLAT), i.e. approximately in the centre of the PV streamer



1471 
shown in Fig.\,\ref{question}. The last parameter (CROT) allows the new grid to be rotated. At the end of this step, the following fields are available on the R20060115\_18 file:\\[0.2cm]



1472 



1473 
\begin{center}



1474 
\begin{tabular}{lll}



1475 
\hline



1476 
Q & Specific humidity ($kg/kg$) & R20060116\_18 \\



1477 
U & (Rotated) Zonal velocity ($m/s$) & R20060116\_18 \\



1478 
V & (Rotated) Meridional velocity ($m/s$) & R20060116\_18 \\



1479 
T & Temperature ($^\circ C$) & R20060116\_18 \\



1480 
P & Pressure ($hPa$) & R20060116\_18 \\



1481 
ORO & Height of topography ($m$) & R20060116\_18 \\



1482 
LAT & Geographical latitude & R20060116\_18 \\



1483 
LON & Geographical longitude & R20060116\_18 \\



1484 
CORIOL & Coriolis parameter ($1/s$) & R20060116\_18 \\



1485 
X & Cartesian distance from centre along x axis ($km$) & R20060116\_18 \\



1486 
Y & Cartesian distance from centre along x axis ($km$) & R20060116\_18 \\



1487 
\hline



1488 
\end{tabular}



1489 
\end{center}



1490 



1491 



1492 



1493 
\subsubsection{Calculating secondary fields [prep3]}



1494 



1495 



1496 
Several additional fields are needed for the PV inversion, the most prominent of course being Ertel's PV itself.



1497 
This field can easily derived from the horizontal wind components, the potential temperature and



1498 
the density. Firstly, the potential temperature $\theta$ is calculated from pressure $p$ (in hPa) and



1499 
temperature $T$ (in K) according to its definition:



1500 



1501 
\[



1502 
\theta = T \cdot (\frac{p_0}{p})^\kappa



1503 
\]



1504 
where the reference pressure is fixed as 1000\,hPa and $\kappa=R/C_p$ is the ratio of the gas constant



1505 
for dry air ($R$) and the specific heat at constant pressure ($c_p$). Numerically, this value



1506 
is in good approximation: $\kappa=0.286$.



1507 
Note that an inherent property of the atmosphere is its stability. This is expressed physically by



1508 
a potential temperature which increases with height. Regions where potential temperature decreases



1509 
with height are hydrostatically unstable and the prerequisites for a PV inversion are not fulfilled



1510 
there. \\



1511 



1512 



1513 
\noindent



1514 
The density $\rho$ (in kg/m$^3$) of the air is determined by the ideal gas equation for air.\\[0.3cm]



1515 



1516 
\[



1517 
\rho = \frac{p}{R_D \cdot T}



1518 
\]



1519 



1520 
\vspace{0.3cm}



1521 
\noindent



1522 
where $R_d$ is the gas constant for dry air and p and T are the pressure (in Pa) and temperature



1523 
(in K), respectively.\\



1524 



1525 
\noindent



1526 
With this definition of potential temperature and density, Ertel's PV is readily calculated according to:



1527 



1528 
\[



1529 
PV = \frac{1}{\rho} \cdot ( \vec{\xi} + f \cdot \hat{k} ) \cdot \vec{\nabla}{\theta}



1530 
\]



1531 



1532 
\vspace{0.3cm}



1533 
\noindent



1534 
Here, $\hat{k}$ is a unit vector pointing in the vertical direction, $f$ is the planetary vorticity



1535 
(typically $10^{4} s^{1}$ in the extratropics),



1536 
i.e. the vorticity due to the Earth's rotation, and $\vec{\xi}$ is the curl $\vec{\nabla} \times \vec{u}$



1537 
of the wind vector $\vec{u}$. In largescale dynamics vertical wind components can be neglected in the



1538 
above formula. The final expression for Ertel's PV then reduces to:



1539 



1540 
\[



1541 
PV = \frac{1}{\rho} \cdot [ ( \frac{\partial v}{\partial x}  \frac{\partial u}{\partial y} + f )



1542 
\cdot \frac{\partial \theta}{\partial z}



1543 
+ \frac{\partial u}{\partial z} \cdot \frac{\partial \theta}{\partial y}



1544 
 \frac{\partial v}{\partial z} \cdot \frac{\partial \theta}{\partial x}]



1545 
\]



1546 



1547 
\vspace{0.3cm}



1548 
\noindent



1549 
In many cases the latter two terms (involving vertical derivatives of horizontal velocity) can



1550 
also be neglected. Here, we keep these terms in order to have a higherorder approximation to Ertel's



1551 
PV. Note that the first term on the right side includes the aforementioned vertical



1552 
derivative of potential temperature. The stability criterion forces this derivative to be positive.



1553 
The theory of atmospheric turbulence leads to a still stronger constraint. Indeed, the whole PV



1554 
must be positive (on the northern hemisphere), since otherwise symmetric instability will set in and



1555 
turbulently adjust the atmosphere to a state where PV is nonnegative.\\



1556 



1557 
\begin{figure}[t]



1558 
\begin{center}



1559 
\begin{minipage}[t]{7.8cm}



1560 
\ForceWidth{1.0\textwidth}



1561 
\BoxedEPSF{zus_pv_k40.eps}



1562 
\end{minipage}



1563 
\hspace{0cm}



1564 
\begin{minipage}[t]{7.8cm}



1565 
\ForceWidth{1.0\textwidth}



1566 
\BoxedEPSF{zus_th_k40.eps}



1567 
\end{minipage}



1568 
\\[0cm]



1569 
\begin{minipage}[t]{7.8cm}



1570 
\ForceWidth{1.0\textwidth}



1571 
\BoxedEPSF{zus_nsq_k40.eps}



1572 
\end{minipage}



1573 
\hspace{0cm}



1574 
\begin{minipage}[t]{7.8cm}



1575 
\ForceWidth{1.0\textwidth}



1576 
\BoxedEPSF{zus_rho_k40.eps}



1577 
\end{minipage}



1578 
\\[0cm]



1579 
\caption{\it Secondary fields:



1580 
Ertel's potential vorticity (upperleft, in $pvu$), potential temperature (upperright, in K),



1581 
squared BruntV\"ais\"al\"a frequency (lowerleft, in $10^{4} s^{2}$) and density (lowerright,



1582 
in $kg m^{3}$). The stratification is a measure of the vertical stratification of the atmosphere.



1583 
High values correspond



1584 
to strong stratification and suppression of vertical motions. Note how the stratospheric PV streamer is



1585 
associated with enhanced stratification. The density, on the other hand, is slightly decreased within the



1586 
streamer's region. All



1587 
fields are plotted on model level 40, corresponding to a geometrical height of approximately 8\,km.



1588 
}



1589 
\label{secfield}



1590 
\end{center}



1591 
\end{figure}



1592 



1593 
\noindent



1594 
The last field which is needed for the PV inversion is strongly related to the vertical derivative



1595 
of potential temperature. It is the socalled squared BruntV\"as\"ala frequency (or stratification),



1596 
defined by



1597 



1598 
\[



1599 
N^2 = \frac{g}{\theta} \cdot \frac{\partial \theta}{\partial z}



1600 
\]



1601 



1602 
\vspace{0.3cm}



1603 
\noindent



1604 
where $g$ is the Earth's gravity. Following the above discussion, a necessary condition for



1605 
atmospheric stability is that $N^2$ is nonnegative. Note that strong stratifications counteract



1606 
vertical motions. Particularly strong stratifications are typically found in the stratosphere,



1607 
($N^2 \approx 10^{4} s^{2}$) whereas the troposphere is considerably weaker stratified



1608 
($N^2 \approx 10^{4} s^{2}$).\\



1609 



1610 
\noindent



1611 
All these secondary fields are calculated with the call to \\[0.3cm]



1612 
\begin{center}



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



1614 
\begin{verbatim}



1615 
inversion.sh prep3



1616 
\end{verbatim}



1617 
\end{minipage}



1618 
\end{center}



1619 



1620 
\vspace{0.3cm}



1621 
\noindent



1622 
After this step, four additional variables are available on the file R20060116\_18. They are listed in the following table and are shown in Fig.\,\ref{secfield}:\\



1623 



1624 



1625 
\begin{center}



1626 
\begin{tabular}{lll}



1627 
\hline



1628 
TH & Potential temperature ($K$) & R20060116\_18 \\



1629 
PV & Ertel's potential vorticity ($pvu$) & R20060116\_18 \\



1630 
NSQ & Squared BruntVai\"al\"a frequency ($1/s^2$) & R20060116\_18 \\



1631 
RHO & Air density ($kg/m^3$) & R20060116\_18 \\



1632 
\hline



1633 
\end{tabular}



1634 
\end{center}



1635 
\vspace{0.5cm}



1636 



1637 
\subsubsection{Identification of a PV anomaly [prep4]}



1638 



1639 



1640 
Having calculated Ertel's PV on height levels, it is now time to specify the modified PV distribution



1641 
which should be obtained by the PV inversion. In short, a PV anomaly has to be identified from the original



1642 
PV field and then the effect of this anomaly on the potential temperature and velocity has to be determined.\\



1643 



1644 
\noindent



1645 
There is no unique and foolproof method how to specify an anomaly. Indeed, this step strongly depends on



1646 
the features which should be studied. In the present example, we look at a socalled stratospheric PV streamer and the aim



1647 
is to remove this streamer from the PV field and thereby to quantify its influence on the stability and wind



1648 
fields in the middle and lower troposphere. Figure\,ref{box} shows the PV in a horizontal crosssection at 8\,km



1649 
height. The PV streamer is clearly discernible as a



1650 
southward intrusion of highPV, i.e. stratospheric air toward the south. In the vertical the streamer reaches



1651 
down to levels of 6\,km (not shown).\\



1652 



1653 
\noindent



1654 
The extraction of the streamer is performed by filtering the original PV field. This filtering is done



1655 
by the call\\[0.3cm]



1656 
\begin{center}



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



1658 
\begin{verbatim}



1659 
inversion.sh prep4



1660 
\end{verbatim}



1661 
\end{minipage}



1662 
\end{center}



1663 



1664 
\vspace{0.3cm}



1665 
\noindent



1666 
Essentially, the program performs



1667 
the following two steps: Choose a grid point within the filtering box (see Fig\,ref{box}), and



1668 
replace the PV value at this grid point by the zonal mean of all PV values. Formally this might be



1669 
written as:



1670 



1671 
\[



1672 
PV \quad \rightarrow \quad <PV>



1673 
\]



1674 



1675 
\vspace{0.3cm}



1676 
\noindent



1677 
where $<>$ denotes the filtering operator which applies only within the filtering box.



1678 
The anomaly could now be determined as the difference of the original and filtered PV:



1679 
It turns out that this simple difference $AN=PV<PV>$ would produce not only the desired positive anomaly associated



1680 
with the streamer, but also surrounding negative PV anomalies. We neglect all these negative anomalies



1681 
by simply setting them to zero. \\[0.2cm]



1682 



1683 
\[



1684 
AN = PV  <PV> \quad \mbox{and} \quad AN \rightarrow AN^+.



1685 
\]



1686 



1687 
\begin{figure}[t]



1688 
\begin{center}



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



1690 
\ForceWidth{1.0\textwidth}



1691 
\BoxedEPSF{box.eps}



1692 
\end{minipage}



1693 
\\



1694 
\caption{\it Ertel's PV on level $k=40$ (corresponding to a height of 8\,km). Additionally the horizontal



1695 
grid is shown, and the filtering box is marked with the bold black lines.}



1696 
\label{box}



1697 
\end{center}



1698 
\end{figure}



1699 



1700 
\vspace{0.3cm}



1701 
\noindent



1702 
The outcome $AN^+$ is a "wellbehaved" PV anomaly, which in turn is used to finally define the modified $PV^M$, i.e.



1703 
the PV field which should be obtained by the inversion:\\



1704 



1705 
\[



1706 
PV^M = PV  AN^+



1707 
\]



1708 



1709 
\vspace{0.3cm}



1710 
\noindent



1711 
The original $PV$, modified $PV^M$ and the anomaly $AN^+$ is shown in Fig.\,\ref{ano} at several model levels. Note how the PV streamer of Fig.\,\ref{question} is well captured by the PV anomaly



1712 



1713 



1714 
\begin{figure}[t]



1715 
\begin{center}



1716 
\begin{minipage}[t]{5cm}



1717 
\ForceWidth{1.0\textwidth}



1718 
\BoxedEPSF{pv_org_k20.eps}



1719 
\end{minipage}



1720 
\hspace{0.0cm}



1721 
\begin{minipage}[t]{5cm}



1722 
\ForceWidth{1.0\textwidth}



1723 
\BoxedEPSF{pv_mod_k20.eps}



1724 
\end{minipage}



1725 
\hspace{0.0cm}



1726 
\begin{minipage}[t]{5cm}



1727 
\ForceWidth{1.0\textwidth}



1728 
\BoxedEPSF{pv_ano_k20.eps}



1729 
\end{minipage}



1730 
\\



1731 
\begin{minipage}[t]{5cm}



1732 
\ForceWidth{1.0\textwidth}



1733 
\BoxedEPSF{pv_org_k30.eps}



1734 
\end{minipage}



1735 
\hspace{0.0cm}



1736 
\begin{minipage}[t]{5cm}



1737 
\ForceWidth{1.0\textwidth}



1738 
\BoxedEPSF{pv_mod_k30.eps}



1739 
\end{minipage}



1740 
\hspace{0.0cm}



1741 
\begin{minipage}[t]{5cm}



1742 
\ForceWidth{1.0\textwidth}



1743 
\BoxedEPSF{pv_ano_k30.eps}



1744 
\end{minipage}



1745 
\\



1746 
\begin{minipage}[t]{5cm}



1747 
\ForceWidth{1.0\textwidth}



1748 
\BoxedEPSF{pv_org_k40.eps}



1749 
\end{minipage}



1750 
\hspace{0.0cm}



1751 
\begin{minipage}[t]{5cm}



1752 
\ForceWidth{1.0\textwidth}



1753 
\BoxedEPSF{pv_mod_k40.eps}



1754 
\end{minipage}



1755 
\hspace{0.0cm}



1756 
\begin{minipage}[t]{5cm}



1757 
\ForceWidth{1.0\textwidth}



1758 
\BoxedEPSF{pv_ano_k40.eps}



1759 
\end{minipage}



1760 
\\



1761 
\begin{minipage}[t]{5cm}



1762 
\ForceWidth{1.0\textwidth}



1763 
\BoxedEPSF{pv_org_k50.eps}



1764 
\end{minipage}



1765 
\hspace{0.0cm}



1766 
\begin{minipage}[t]{5cm}



1767 
\ForceWidth{1.0\textwidth}



1768 
\BoxedEPSF{pv_mod_k50.eps}



1769 
\end{minipage}



1770 
\hspace{0.0cm}



1771 
\begin{minipage}[t]{5cm}



1772 
\ForceWidth{1.0\textwidth}



1773 
\BoxedEPSF{pv_ano_k50.eps}



1774 
\end{minipage}



1775 
\\



1776 
\caption{\it Original (left), modified (middle) and anomaly (right) Ertel's PV for 16 January 2006, 18\,UTC. The different rows correspond to the levels (from top to bottom) 20, 30 40 50 corresponding to the heights 4, 6, 8 and 10\,km.}



1777 
\label{ano}



1778 
\end{center}



1779 
\end{figure}



1780 



1781 
\noindent



1782 
The filtering



1783 
is characterised by the values in six lines of the parameter file. \\[0.3cm]



1784 
\begin{center}



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



1786 
\begin{verbatim}



1787 
BEGIN ANOMALY



1788 
BOX_XMIN = 1200.;



1789 
BOX_XMAX = 1200.;



1790 
BOX_YMIN = 1200.;



1791 
BOX_YMAX = 1200.;



1792 
BOX_ZMIN = 2000.;



1793 
BOX_ZMAX = 10000.;



1794 
NFILTER = 5 ;



1795 
BOUND_XY = 500.;



1796 
BOUND_Z = 500.;



1797 
END ANOMALY



1798 
\end{verbatim}



1799 
\end{minipage}



1800 
\end{center}



1801 



1802 
\vspace{0.6cm}



1803 
\noindent



1804 
The lines specify the filtering box in the westeast direction (from 1200\,km to



1805 
1200\,km, from BOX\_XMIN to BOX\_XMAX), in the southnorth direction (from 1200\,km to 1200\,km, from BOX\_YMIN to BOX\_YMAX), and in the vertical



1806 
(from 2000\,m to 10000\,m, from BOX\_ZMIN to BOX\_ZMAX). The filtering is not only applied once, but



1807 
(as given by NFILTER) $n=5$ times.



1808 
Finally the last two parameters describe the handling at the boundaries of the filtering box. In fact,



1809 
it is advantageous to force the anomaly PV to zero at the boundaries of the box. This is accomplished



1810 
by multiplying the PV anomaly with a function f$_{\partial F}$ which is 0 without the box, then monotonically increases



1811 
to 1 within a boundary zone, and then is constant 1 in the interior of the box. The width of the boundary



1812 
zone is given as 500\,km in the horizontal direction (BOUND\_XY) and 500\,m in the vertical direction (BOUND\_Z).\\



1813 



1814 
\noindent



1815 
The modified PV field and the anomaly are written to the R20060116\_18 file.



1816 
Additionally the boundary values are written. The velocity fields must be specified on the lateral



1817 
boundaries, whereas the potential temperature anomaly must be specified on the lower and upper boundary. In



1818 
the present settings, all boundary values are zero. The following table lists all fields which are additionally



1819 
written to the R file.



1820 



1821 
\begin{center}



1822 
\begin{tabular}{lll}



1823 
\hline



1824 
PV\_FILT & Modified Ertel's PV ($pvu$) & R20060116\_18 \\



1825 
PV\_ANOM & Anomaly in Ertel' PV ($pvu$) & R20060116\_18 \\



1826 
UU\_ANOM & Lateral boundary condition of zonal wind($m/s$) & R20060116\_18 \\



1827 
VV\_ANOM & Lateral boundary condition of meridional wind ($m/s$) & R20060116\_18 \\



1828 
TH\_ANOM & Upper and lower boundary of potential temperature ($K$) & R20060116\_18 \\



1829 
\hline



1830 
\end{tabular}



1831 
\end{center}



1832 
\vspace{0.5cm}



1833 



1834 



1835 
\subsubsection{Splitting into different files [prep5]}



1836 



1837 
\begin{figure}[t]



1838 
\begin{center}



1839 
\begin{minipage}[t]{11.5cm}



1840 
\hspace{2.5cm}



1841 
\ForceWidth{1.0\textwidth}



1842 
\BoxedEPSF{splitting.eps}



1843 
\end{minipage}



1844 
\\



1845 
\caption{\it Splitting of the initial R20060116\_18 file into four different files with



1846 
prefix ORG, MOD, ANO and REF. }



1847 
\label{split}



1848 
\end{center}



1849 
\end{figure}



1850 



1851 
\noindent



1852 
So far, all new fields are written to the R20060116\_18 file. In the further steps it is



1853 
advantageous to split this file into several files. One file should get all original



1854 
fields, one only the modified fields and one only anomaly fields. Additionally, some



1855 
"static" fields, such as orography and Coriolis parameter, should be written to a special



1856 
file, from here on called the reference file.\\



1857 



1858 
\noindent



1859 
The splitting is done with the call\\[0.3cm]



1860 
\begin{center}



1861 
\begin{minipage}[t]{15cm}



1862 
\begin{verbatim}



1863 
inversion.sh prep4



1864 
\end{verbatim}



1865 
\end{minipage}



1866 
\end{center}



1867 



1868 
\vspace{0.3cm}



1869 
\noindent



1870 
The result of this step is a set of new netcdf files with prefix ORG, MOD, ANO and REF. The following table



1871 
(and Fig.\,\ref{split}) illustrates which variables are written to which file. Note that some variables are renamed in this process.



1872 
\\



1873 



1874 
\begin{center}



1875 
\begin{tabular}{rrrr}



1876 
\hline



1877 
Old variable name & New variable name & Source file & Destination file \\



1878 
\hline



1879 
PV & PV & R20060115\_18 & ORG\_20060116\_18 \\



1880 
U & U & R20060115\_18 & ORG\_20060116\_18 \\



1881 
V & V & R20060115\_18 & ORG\_20060116\_18 \\



1882 
TH & TH & R20060115\_18 & ORG\_20060116\_18 \\



1883 
Q & Q & R20060115\_18 & ORG\_20060116\_18 \\



1884 
P & P & R20060115\_18 & ORG\_20060116\_18 \\



1885 
T & T & R20060115\_18 & ORG\_20060116\_18 \\



1886 
\hline



1887 
PV\_FILT & PV\_AIM & R20060115\_18 & MOD\_20060116\_18 \\



1888 
U & U & R20060115\_18 & MOD\_20060116\_18 \\



1889 
V & V & R20060115\_18 & MOD\_20060116\_18 \\



1890 
Q & Q & R20060115\_18 & MOD\_20060116\_18 \\



1891 
P & P & R20060115\_18 & MOD\_20060116\_18 \\



1892 
T & T & R20060115\_18 & MOD\_20060116\_18 \\



1893 
NSQ & NSQ & R20060115\_18 & MOD\_20060116\_18 \\



1894 
RHO & RHO & R20060115\_18 & MOD\_20060116\_18 \\



1895 
TH & TH & R20060115\_18 & MOD\_20060116\_18 \\



1896 
\hline



1897 
\end{tabular}



1898 



1899 
\begin{tabular}{rrrr}



1900 
\hline



1901 
Old variable name & New variable name & Source file & Destination file \\



1902 
\hline



1903 
PV\_ANOM & PV & R20060115\_18 & ANO\_20060115\_18 \\



1904 
TH\_ANOM & TH & R20060115\_18 & ANO\_20060115\_18 \\



1905 
UU\_ANOM & U & R20060115\_18 & ANO\_20060115\_18 \\



1906 
VV\_ANOM & V & R20060115\_18 & ANO\_20060115\_18 \\



1907 
\hline



1908 
ORO & ORO & R20060115\_18 & REF\_20060115\_18 \\



1909 
X & X & R20060115\_18 & REF\_20060115\_18 \\



1910 
Y & Y & R20060115\_18 & REF\_20060115\_18 \\



1911 
LAT & LAT & R20060115\_18 & REF\_20060115\_18 \\



1912 
LON & LON & R20060115\_18 & REF\_20060115\_18 \\



1913 
CORIOL & CORIOL & R20060115\_18 & REF\_20060115\_18 \\



1914 
\hline



1915 
\end{tabular}



1916 
\end{center}



1917 
\vspace{0.5cm}



1918 



1919 



1920 



1921 
\subsubsection{Specifying the reference profile [prep6]}



1922 



1923 



1924 
A reference vertical profile must be



1925 
defined for the PV inversion. The profile should be as near as possible to the real profile. Any quantity is defined in



1926 
the inversion as a deviation from this profile. Hence if the profile is well chosen only small



1927 
deviations can be expected and the assumed linearisation (see later) is better. Physically we determine



1928 
a good reference profile by taking the area mean over the subdomain, i.e at every model level the



1929 
mean over potential temperature $\theta$, squared Brunt Vais\"al\"a frequency $N^2$ and density $\rho$



1930 
is calculated. The resulting profiles are shown in Fig.\,\ref{reference} as a function of height, and additionally



1931 
as a function of pressure.\\



1932 



1933 
\noindent



1934 
Note how the density of the reference profile exponentially decreases



1935 
decreases with increasing height. If pressure is used as vertical coordinate, the decrease in density



1936 
is essentially linear with decreasing pressure. The profile of potential temperature increases with



1937 
increasing height: This is a necessary condition for the convective stability of the air column. In fact,



1938 
any region where potential temperature decreases with increasing height is unstable and characterised by high levels



1939 
of turbulence. The stability expressed by this potential temperature increase with height is numerically



1940 
determined by the socalled (squared) BruntVais\"al\"a frequency, shown in the left panels. Since this



1941 
field directly goes into the solution of the quasigeostrophic PV equation (see section~3.5), it must



1942 
be guaranteed that it always remains positive.\\



1943 



1944 
\begin{figure}[t]



1945 
\begin{center}



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



1947 
\ForceWidth{1.0\textwidth}



1948 
\BoxedEPSF{reference.eps}



1949 
\end{minipage}



1950 
\\



1951 
\caption{\it Reference profiles of potential temperature $\theta$, squared Brunt Vais\"al\"a frequency $N^2$ and density $\rho$ as a function of height (upper panels) and pressure (lower panels). The



1952 
reference profiles correspond to the areaaveraged values over the inversion subdomain. }



1953 
\label{reference}



1954 
\end{center}



1955 
\end{figure}



1956 



1957 
\noindent



1958 
The reference profile is calculated with the call to \\[0.3cm]



1959 
\begin{center}



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



1961 
\begin{verbatim}



1962 
inversion.sh prep5



1963 
\end{verbatim}



1964 
\end{minipage}



1965 
\end{center}



1966 



1967 
\vspace{0.3cm}



1968 
\noindent



1969 
After his step some additional variables are available on the REF file. These are:\\



1970 



1971 
\begin{center}



1972 
\begin{tabular}{lll}



1973 
\hline



1974 
NSQREF & Reference squared BruntVais\"al\"a frequency ($1/s^2$) & MOD\_20060116\_18 \\



1975 
RHOREF & Reference air density ($kg/m^3$) & MOD\_20060116\_18 \\



1976 
PREREF & Reference pressure ($hPa$) & MOD\_20060115\_18 \\



1977 
THETAREF & Reference potential temperature ($K$) & MOD\_20060116\_18 \\



1978 
ZREF & Reference geopotential height ($m$) & MOD\_20060116\_18 \\



1979 
\hline



1980 
\end{tabular}



1981 
\end{center}



1982 
\vspace{0.5cm}



1983 



1984 
\subsubsection{Adding the coastlines [prep7]}



1985 



1986 
\begin{figure}[t]



1987 
\begin{center}



1988 
\begin{minipage}[t]{7.8cm}



1989 
\ForceWidth{1.0\textwidth}



1990 
\BoxedEPSF{coast1.eps}



1991 
\end{minipage}



1992 
\hspace{0.0cm}



1993 
\begin{minipage}[t]{7.8cm}



1994 
\ForceWidth{1.0\textwidth}



1995 
\BoxedEPSF{coast2.eps}



1996 
\end{minipage}



1997 
\\



1998 
\caption{\it Coastlines transformed into the quasicartesian coordinate system. In the left panel, the



1999 
coastlines are shown in a system with rotated longitude and latitude as axes, in the right panel they



2000 
are shown in a x/y coordinate system.}



2001 
\label{coastline}



2002 
\end{center}



2003 
\end{figure}



2004 



2005 
At this stage it is not possible to plot the continental coastlines in the quasicartesian coordinate frame.



2006 
In fact, this becomes only feasible if the geographical coastlines (as given in latitude/longitude coordinates) are transformed into the new quasicartesian coordinate system. This is done



2007 
with the call\\[0.3cm]



2008 
\begin{center}



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



2010 
\begin{verbatim}



2011 
inversion.sh prep7



2012 
\end{verbatim}



2013 
\end{minipage}



2014 
\end{center}



2015 



2016 
\vspace{0.3cm}



2017 
\noindent



2018 
In Fig.\,\ref{coastline} the transformed coordinates are shown. On the left panel, the coordinate axis are the rotated longitude and latitude and on the left panel they are given by the distance (in km) from the coordinate's center (see Fig.\,14).



2019 
In this step the following additional fields are written to the REF\_20060116\_18 file:



2020 



2021 
\begin{center}



2022 
\begin{tabular}{lll}



2023 
\hline



2024 
COAST\_LON & Longitude of coastline & REF\_20060115\_18 \\



2025 
COAST\_LAT & Latitude of coastline & REF\_20060115\_18 \\



2026 
COAST\_RLON & Rotated longitude of coastline & REF\_20060115\_18 \\



2027 
COAST\_RLAT & Rotated latitude of coastline & REF\_20060115\_18 \\



2028 
COAST\_X & X coordinate of coastline & REF\_20060115\_18 \\



2029 
COAST\_Y & Y coordinate of coastline & REF\_20060115\_18 \\



2030 
\hline



2031 
\end{tabular}



2032 
\end{center}



2033 



2034 
\vspace{0.3cm}



2035 
\noindent



2036 
The coastline in geographical coordinates is taken from a global data set which can be retrieved from



2037 
the National Geophysical Data Center



2038 
(rimmer.ngdc.noaa.gov/coast/). Only, the part of the coastlines inside the rotated coordinate system are specified. This allows



2039 
an efficient plotting of the coastlines. If Matlab is used as a visualisation tool, the plotting can be



2040 
done with the following few lines, depending whether rotated longitude/latitude or X/Y is used as coordinate



2041 
axes.\\[0.0cm]



2042 
\begin{center}



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



2044 
\begin{verbatim}



2045 
%Read input file



2046 
inp = ncget('REF_20060116_18');



2047 



2048 
% Plot continental boundaries (rotated latitude/longitude)



2049 
figure;



2050 
plot(inp.COAST_RLON.data,inp.COAST_RLAT.data,'k');



2051 



2052 
% Plot continental boundaries (x/y)



2053 
figure;



2054 
plot(inp.COAST_X.data,inp.COAST_Y.data,'k');



2055 
\end{verbatim}



2056 
\end{minipage}



2057 
\end{center}



2058 



2059 
\vspace{0.6cm}



2060 
\noindent



2061 
The first command reads the reference file. All variables and parameters are then available in the Matlab structure inp. The COAST subfields of this structure can then immediately be plotted. Sometimes it turns



2062 
out to be better to set all coordinates outside the plotting domain to NaN.



2063 



2064 



2065 
\subsubsection{Moving the files to the run directory [prep8]}



2066 



2067 
At this stage, all essential preparatory steps for the PV inversion are done: The modified PV field is defined, the



2068 
boundary conditions are written and the reference profile is available. In this last step, the files are moved from the



2069 
input directory to the run directory, where all iterations of the PV inversions are performed. The moving of the



2070 
files is done with\\[0.3cm]



2071 
\begin{center}



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



2073 
\begin{verbatim}



2074 
inversion.sh prep8



2075 
\end{verbatim}



2076 
\end{minipage}



2077 
\end{center}



2078 



2079 
\vspace{0.3cm}



2080 
\noindent



2081 
The ORG, MOD, ANO and REF files are now available in the run directory (see Fig.\,8). No further processes are



2082 
run in the input directory.



2083 



2084 



2085 



2086 
\subsection{Inversion  Iterative steps toward modified PV field [pvin]}



2087 



2088 
\begin{figure}[t]



2089 
\begin{center}



2090 
\begin{minipage}[t]{14cm}



2091 
\ForceWidth{1.0\textwidth}



2092 
\BoxedEPSF{iterations.eps}



2093 
\end{minipage}



2094 
\\



2095 
\caption{\it Flowchart of the different inversion steps. Additionally, to the right of the different blocks, the



2096 
"flow" of information is illustrated. For instance, MOD $\rightarrow$ ANO means that the input information



2097 
is taken from the MOD file and then written to the ANO file. }



2098 
\label{iterations}



2099 
\end{center}



2100 
\end{figure}



2101 



2102 
The following steps build the core of the PV inversion. An interior anomaly of Ertel's PV



2103 
and the boundary values for the inversion were specified in the previous preparatory steps.



2104 
Now the flow anomalies (wind vectors, temperature, pressure and potential temperature) associated



2105 
with this PV anomaly has to be determined. This is the aim of the inversion. More specifically



2106 
the process must be split again into several distinct steps. These different steps are shown as a flowchart in



2107 
Fig.\,\ref{iterations}.



2108 



2109 



2110 



2111 



2112 
\subsubsection{Calculation of secondary fields [pvin1]}



2113 



2114 
Secondary fields have to be calculated. These fields are potential temperature, potential vorticity



2115 
density and squared BruntVais\"al\"a frequency. The details of these calculations were already presented



2116 
in section~4.4.3, and are not repeated here. All the secondary fields are calculated with the call\\[0.3cm]



2117 
\begin{center}



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



2119 
\begin{verbatim}



2120 
inversion.sh pvin1



2121 
\end{verbatim}



2122 
\end{minipage}



2123 
\end{center}



2124 



2125 
\vspace{0.3cm}



2126 
\noindent



2127 
The call subsequently calculates potential temperature, squared BruntV\"as\"ala frequency, density and Ertel's PV. Where the height levels are below topography,



2128 
missing data values are written. The following table lists the new fields which are written to the MOD file:



2129 



2130 
\begin{center}



2131 
\begin{tabular}{lll}



2132 
\hline



2133 
TH & Potential temperature ($K$) & MOD\_20060116\_18 \\



2134 
PV & Ertel's potential vorticity ($pvu$) & MOD\_20060116\_18 \\



2135 
NSQ & Squared BruntVai\"al\"a frequency ($1/s^2$) & MOD\_20060116\_18 \\



2136 
RHO & Air density ($kg/m^3$) & MOD\_20060116\_18 \\



2137 
\hline



2138 
\end{tabular}



2139 
\end{center}



2140 
\vspace{0.5cm}



2141 



2142 
\subsubsection{Transforming Ertel's PV to quasigeostrophic PV [pvin2]}



2143 



2144 
\begin{figure}[t]



2145 
\begin{center}



2146 
\begin{minipage}[t]{14.0cm}



2147 
\hspace{1.0cm}



2148 
\ForceWidth{1.0\textwidth}



2149 
\BoxedEPSF{pv_to_qgpv.eps}



2150 
\end{minipage}



2151 
\\



2152 
\caption{\it Determination of Ertel's PV anomaly and its transformation into



2153 
an anomaly of quasigeostrophic PV. The calculated Ertel's PV is subtracted from



2154 
Ertel's PV which should be reached, and then this difference is transformed into



2155 
an approximate quasigeostrophic PV. The input is taken from the MOD file, and the



2156 
output is written to the ANO file.}



2157 
\label{pvtoqg}



2158 
\end{center}



2159 
\end{figure}



2160 



2161 
So far only Ertel's PV was used. On the other hand, the inverson problem is formulated for the



2162 
quasigeostrophic PV. Therefore a conversion between the two has to be done. In a first approximation



2163 
Ertels's PV is given by (for clarity it is denoted as EPV in this section):\\[0.2cm]



2164 



2165 
\[



2166 
EPV = \frac{1}{\rho} \cdot ( \xi + f ) \cdot \frac{\partial \theta}{\partial z}



2167 
\]



2168 



2169 
\vspace{0.3cm}



2170 
\noindent



2171 
where $\rho$ is the density, $\xi$ the relative vorticity, $f$ the Coriolis parameter and $\theta$ the potential temperature.



2172 
Ertel's PV can be expressed in terms of the reference profile and deviation thereof. If we set



2173 
$\theta=\theta_0 + \theta^*$ and $\rho=\rho_0 + \rho^*$, consider the definition



2174 
of the reference stratification $N^2_0=g/\theta_0 \cdot \partial {\theta_0}/{\partial z}$ and furthermore



2175 
apply the quasigeostrophic assumption $\xi \ll f$, the following approximate



2176 
expression hold in first order (we have neglected the density perturbation $\rho^*$ in the denominator):\\[2mm]



2177 



2178 
\[



2179 
EPV = \frac{\theta_0 N_0^2}{g \rho_0} \cdot (\xi + \frac{f g}{\theta_0 N_0^2} \cdot \frac{\partial \theta^*}{\partial z})



2180 
\]



2181 



2182 
\vspace{0.3cm}



2183 
\noindent



2184 
On the other hand, the quasigeostrophic PV (denoted by qgPV in this section) is defined by the streamfunction $\psi$\\[0.2cm]



2185 



2186 
\[



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



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



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



2190 
\]



2191 



2192 
\vspace{0.3cm}



2193 
\noindent



2194 
where $\rho_0$ and $N_0^2$ denotes a reference profile of density and stratification (as defined in the



2195 
preparatory steps and written to the REF file). The following relationship between the streamfunction and the meteorological fields persists:\\[2mm]



2196 



2197 
\[



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



2199 
\quad \quad



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



2201 
\quad \quad



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



2203 
\]



2204 



2205 
\vspace{0.3cm}



2206 
\noindent



2207 
If these expressions are inserted into the above equation for $q$, one gets:\\[2mm]



2208 



2209 
\[



2210 
qgPV = \xi + \frac{f^2}{\rho_0} \cdot \frac{\partial}{\partial z}



2211 
(\frac{\rho_0 g}{N_0^2 f \theta_0} \cdot \theta^* )



2212 
\]



2213 



2214 
\vspace{0.3cm}



2215 
\noindent



2216 
In a first approximation, we treat the reference profile and the Coriolis parameter as constant. Then,



2217 
this last expression is further simplified to:\\[2mm]



2218 



2219 
\[



2220 
qgPV = \xi + \frac{f g}{N_0^2 \theta_0} \cdot \frac{\partial \theta^*}{\partial z}



2221 
\]



2222 



2223 



2224 
\vspace{0.3cm}



2225 
\noindent



2226 
With this form, it is easy to see that the following approximate relationship between Ertel's and the



2227 
quasigeostrophic PV holds. \\[2mm]



2228 



2229 
\[



2230 
qgPV \approx \frac{\rho_0 g}{\theta_0 N_0^2} \cdot EPV  f



2231 
\]



2232 



2233 
\vspace{0.3cm}



2234 
\noindent



2235 
Note that according to the previous derivation, this equation does not hold exactly. It is based upon



2236 
an assumption of linearisation and on approximate forms of Ertel's PV. Nevertheless, we can expect it



2237 
to be a good firstorder approximation. In the PV inversion tool, use is made of this relationship,



2238 
but since it is only approximate, an iterative technique has to be applied. Finally note that this



2239 
formula can be easily applied to PV anomalies. Indeed, if $\Delta PV$ refers to a anomaly in Ertel's



2240 
PV, the corresponding anomaly in quasigeostrophic PV is given by:\\[2mm]



2241 



2242 
\[



2243 
\Delta(qgPV) \approx \frac{\rho_0 g}{\theta_0 N_0^2} \cdot \Delta(EPV)



2244 
\]



2245 



2246 



2247 
\vspace{0.3cm}



2248 
\noindent



2249 
It is in this latter form that the program call\\[0.3cm]



2250 
\begin{center}



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



2252 
\begin{verbatim}



2253 
inversion.sh pvin2



2254 
\end{verbatim}



2255 
\end{minipage}



2256 
\end{center}



2257 



2258 
\begin{figure}[t]



2259 
\begin{center}



2260 
\begin{minipage}[t]{8.1cm}



2261 
\ForceWidth{1.0\textwidth}



2262 
\BoxedEPSF{qgpv_ano_ew.eps}



2263 
\end{minipage}



2264 
\hspace{0.5cm}



2265 
\begin{minipage}[t]{8.1cm}



2266 
\ForceWidth{1.0\textwidth}



2267 
\BoxedEPSF{qgpv_ano_ns.eps}



2268 
\end{minipage}



2269 
\\



2270 
\caption{\it Quasigeostrophic PV anomaly in an east/west (left) and an north/south (right)



2271 
vertical crosssection for 16 January 2006, 18\,UTC. }



2272 
\label{qgano}



2273 
\end{center}



2274 
\end{figure}



2275 



2276 
\vspace{0.3cm}



2277 
\noindent



2278 
makes the conversion. It takes the ErtelPV from the MOD file, and subtracts it from the PV\_AIM which is also



2279 
available on the MOD file. The difference between the two PV is the new anomaly in Ertel's PV which has to



2280 
be inverted in this iterative step. The approximate relationship between quasigeostrophic PV and Ertel's



2281 
PV is then used to get a corresponding anomaly of the former one (see Fig.\,ref{pvtoqg}).



2282 
Note that this program handles also missing data values. If at some grid point Ertel's PV is not defined,



2283 
i.e. the missing data flag is set, the corresponding value of quasigeostrophic PV is set to zero. It is not



2284 
set to the missing data value since no check will be done by the PV inversion program. Therefore, it is



2285 
best to treat every missing value as a vanishing anomaly of quasigeostrophic PV. A vertical crosssection



2286 
of the anomaly in quasigeostrophic PV is shown in Fig.\,\ref{qgano}.



2287 
The additional or changed fields on the file ANO\_20060116\_18 are given in the following table:\\



2288 



2289 
\begin{center}



2290 
\begin{tabular}{lll}



2291 
\hline



2292 
PV & Ertel's PV anomaly ($pvu$) & ANO\_20060116\_18 \\



2293 
QGPV & Anomaly in quasigeostrophic PV ($s^{1}$) & ANO\_20060116\_18 \\



2294 
\hline



2295 
\end{tabular}



2296 
\end{center}



2297 



2298 



2299 



2300 



2301 



2302 
\subsubsection{Inversion of the quasigeostrophic PV [pvin3]}



2303 



2304 
At this stage the interior distribution of quasigeostrophic PV and the boundary values



2305 
are given, and therefore the inversion can be performed in order to get the streamfunction.



2306 
Numerical aspects of the inversion were already presented in section~2. Here, we focus on the



2307 
practical aspects. The numerical problem solved in this step is expressed in the following



2308 
boundary value problem: \\[0.3cm]



2309 



2310 
\[



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



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



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



2314 
\]



2315 



2316 
\vspace{0.3cm}



2317 
\noindent



2318 
with the predescribed values\\[0.3cm]



2319 



2320 
\[



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



2322 
\quad \quad



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



2324 
\quad \quad



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



2326 
\]



2327 



2328 
\vspace{0.3cm}



2329 
\noindent



2330 
on the upper and lower boundaries (for $\theta^*$) and on the lateral boundaries (for $u$ and $v$). The inversion program is called by\\[0.3cm]



2331 
\begin{center}



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



2333 
\begin{verbatim}



2334 
inversion.sh pvin3



2335 
\end{verbatim}



2336 
\end{minipage}



2337 
\end{center}



2338 



2339 
\vspace{0.3cm}



2340 
\noindent



2341 
Note that the inversion affects only the file ANO\_20060116\_18 because the inversion is



2342 
performed for an anomaly of quasigeostrophic PV.



2343 
The final aim of this step is to calculate the wind, temperature, pressure and potential temperature



2344 
perturbations which are associated with the specified anomaly of quasigeostrophic PV.



2345 
The call to the program yields some information about the ongoing iterative steps. First of all,



2346 
the spectral width in the x, y and zdirection is written to screen. In the present example this



2347 
is:\\[0.3cm]



2348 
\begin{center}



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



2350 
\begin{small}



2351 
\begin{verbatim}



2352 
Spectrum of the matrix in each direction



2353 
Spectrum = 1.028528 0.9136161 1.231250



2354 
\end{verbatim}



2355 
\end{small}



2356 
\end{minipage}



2357 
\end{center}



2358 



2359 
\vspace{0.3cm}



2360 
\noindent



2361 
The spectral widths $\lambda_x,\lambda_y, \lambda_z$ (in the same order as in the above program output) are important to decide whether the inversion can be reasonably be performed. The



2362 
criterion is based upon the matrix elements defining the inversion problem. For a reasonable setting



2363 
the following inequalities should be fulfilled:



2364 



2365 
\[



2366 
\lambda_x > 1/2 \cdot \lambda_y \quad \quad \mbox{and} \quad \quad \lambda_x > 1/2 \cdot \lambda_z



2367 
\]



2368 
\[



2369 
\lambda_y > 1/2 \cdot \lambda_x \quad \quad \mbox{and} \quad \quad \lambda_y > 1/2 \cdot \lambda_z



2370 
\]



2371 
\[



2372 
\lambda_z > 1/2 \cdot \lambda_x \quad \quad \mbox{and} \quad \quad \lambda_z > 1/2 \cdot \lambda_y



2373 
\]



2374 



2375 
\vspace{0.3cm}



2376 
\noindent



2377 
which is evidently fulfilled in our case study (otherwise the program would abort with the error



2378 
message that the grid dimensions are not large enough). Then some runtime statistics is provided for the



2379 
iterative solution of the elliptical partial differential equation. As described in section~2 the



2380 
solution of the quasigeostrophic PV equation is found by means of an successive overrelaxation (SOR)



2381 
technique. Starting from an initial estimate of the streamfunction $\psi$ several hundred iterations of



2382 



2383 



2384 
\[



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



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



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



2388 
\]



2389 



2390 
\vspace{0.3cm}



2391 
\noindent



2392 
are performed (see section~2 for details). Essentially two statistical measures are provided: The first one $\psi_{gauge}$



2393 
measures the amplitude of the streamfunction which is produced in the course of the iteration,



2394 
the second one is a direct check whether the iteration solution converges. In fact, after several iterations the sofar



2395 
reached streamfunction is used to calculate the quasigeostrophic $PV^{c}$ according to the formula:\\[0.3cm]



2396 



2397 
\[



2398 
PV^c = \frac{\partial^2\psi}{\partial x^2} + \frac{\partial^2\psi}{\partial y^2} +



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



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



2401 
\]



2402 



2403 
\vspace{0.3cm}



2404 
\noindent



2405 
This value can then be compared to the quasigeostrophic



2406 
$PV^{i}$ distribution which was given as input to the inversion algorithm. The following norm $\Delta^2$ is then



2407 
calculated as a measure of the convergence:\\[0.3cm]



2408 



2409 
\[



2410 
\Delta^2 = \sum_{i,j,k} (\psi^{i}_{i,j,k}\psi^{c}_{i,j,k})^2



2411 
\]



2412 



2413 
\vspace{0.3cm}



2414 
\noindent



2415 
Some sample output is reproduced below for the case study, the first column corresponding to $\psi_{gauge}$ and the



2416 
second one to $\Delta^2$:\\[0.3cm]



2417 
\begin{center}



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



2419 
\begin{verbatim}



2420 
psigauge = 0.000E+00 deltasq = 0.850E09



2421 
psigauge = 0.000E+00 deltasq = 0.112E09



2422 
psigauge = 0.681E+03 deltasq = 0.643E10



2423 
psigauge = 0.681E+03 deltasq = 0.466E10



2424 
psigauge = 0.806E+03 deltasq = 0.367E10



2425 
psigauge = 0.806E+03 deltasq = 0.303E10



2426 
psigauge = 0.547E+03 deltasq = 0.259E10



2427 
psigauge = 0.547E+03 deltasq = 0.228E10



2428 
psigauge = 0.398E+03 deltasq = 0.206E10



2429 
psigauge = 0.398E+03 deltasq = 0.191E10



2430 
psigauge = 0.309E+03 deltasq = 0.180E10



2431 
\end{verbatim}



2432 
\end{minipage}



2433 
\end{center}



2434 



2435 
\vspace{0.6cm}



2436 
\noindent



2437 
Note that the gauge streamfunction reaches a value of order 3000 (corresponding to a substantial anomaly of wind



2438 
and temperature), and that the convergence measure $\Delta^2$ decreases



2439 
by an order of magnitude from the beginning to the end. At the end of this step, the streamfunction associated with the quasigeostrophic PV anomaly is determined. It is shown in Fig.\,\ref{stream} at three horizontal crosssections. The



2440 
minimum of the streamfunction is found in the interior of the domain, in agreement with the cyclonic flow around the



2441 
positive anomaly of quasigeostrophic PV. \\



2442 



2443 
\begin{figure}[t]



2444 
\begin{center}



2445 
\begin{minipage}[t]{5cm}



2446 
\ForceWidth{1.0\textwidth}



2447 
\BoxedEPSF{psi_ano_k01.eps}



2448 
\end{minipage}



2449 
\hspace{0.0cm}



2450 
\begin{minipage}[t]{5cm}



2451 
\ForceWidth{1.0\textwidth}



2452 
\BoxedEPSF{psi_ano_k20.eps}



2453 
\end{minipage}



2454 
\hspace{0.0cm}



2455 
\begin{minipage}[t]{5cm}



2456 
\ForceWidth{1.0\textwidth}



2457 
\BoxedEPSF{psi_ano_k40.eps}



2458 
\end{minipage}



2459 
\\



2460 
\caption{\it Streamfunction associated with the quasigeostrophic PV anomaly shown in Fig.\,ref{qgano}. The three panels show horizontal crosssections at 0, 4\,km and 8\,km above ground.}



2461 
\label{stream}



2462 
\end{center}



2463 
\end{figure}



2464 



2465 
\noindent



2466 
Now, additional perturbation fields can be computed from this streamfunction: Potential temperature,



2467 
wind components in x and y direction, pressure, and temperature. The corresponding formula are:\\[0.5cm]



2468 



2469 
\begin{eqnarray*}



2470 
\theta^* & = & \frac{f}{g} \cdot \theta_0 \cdot \frac{\partial \psi}{\partial z}\\[0.2cm]



2471 
u & = &  \frac{\partial \psi}{\partial y}\\



2472 
\end{eqnarray*}



2473 
\begin{eqnarray*}



2474 
v & = & \frac{\partial \psi}{\partial x}\\



2475 
\end{eqnarray*}



2476 



2477 
\vspace{0.3cm}



2478 
\noindent



2479 
The anomaly of potential temperature $\theta^*$ and the two horizontal wind components $u$ and $v$ can directly determined from the



2480 
streamfunction. Moreover, the pressure anomaly is directly proportional to the streamfunction, which is only



2481 
determined up to an additive constant. In order to obtain a vanishing pressure at the boundary of the



2482 
inversion domain, we firstly calculate such a constant $\psi_0$ (essentially the average over the domain boundary) and then set the pressure to:\\[0.2cm]



2483 



2484 
\[



2485 
p^* = \rho_0 \cdot f \cdot (\psi  \psi_0)



2486 
\]



2487 



2488 
\vspace{0.3cm}



2489 
\noindent



2490 
The temperature anomaly is now easily calculated from pressure and potential temperature according the relationship



2491 
$\theta = T \cdot (p_ref/p)^\kappa$, where $T$ is the temperature in Kelvin, $p_ref$ is the reference pressure level (1000\,hPa) and $\kappa$ is the ratio between the gas constant for dry air ($R_d$) and the specific heat at constant pressure ($c_p$). Consistent with previous approximations we set:\\[0.2cm]



2492 



2493 
\[



2494 
T^* = (\frac{p_0}{p_ref})^\kappa \cdot (\theta^* + \kappa \cdot \theta_0 \cdot \frac{p^*}{p_0} )



2495 
\]



2496 



2497 
\vspace{0.3cm}



2498 
\noindent



2499 
Figure\,\ref{anofield} shows some of these perturbation fields. The positive anomaly of quasigeostrophic PV is



2500 
associated with a cyclonic (anticlockwise) flow field, which reaches down to surface levels. Furthermore, we know from theoretical studies that a positive PV anomaly goes along with an increase of potential temperature above the anomaly, and correspondingly a decrease of



2501 
potential temperature below. This is nicely confirmed in the perturbation field of the calculation.



2502 
Finally, the cyclonic flow must be associated with a local pressure minimum, according to the



2503 
geostrophic wind equation. This is also fulfilled in our calculation. The following table gives all anomaly fields on the ANO\_10060116\_18 file which are calculated:



2504 



2505 
\begin{figure}[t]



2506 
\begin{center}



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



2508 
\ForceWidth{1.0\textwidth}



2509 
\BoxedEPSF{ano_00_k35_pr.eps}



2510 
\end{minipage}



2511 
\hspace{0.6cm}



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



2513 
\ForceWidth{1.0\textwidth}



2514 
\BoxedEPSF{ano_00_k35_tt.eps}



2515 
\end{minipage}



2516 
\\[0.3cm]



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



2518 
\ForceWidth{1.0\textwidth}



2519 
\BoxedEPSF{ano_00_k35_uv.eps}



2520 
\end{minipage}



2521 
\hspace{0.6cm}



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



2523 
\ForceWidth{1.0\textwidth}



2524 
\BoxedEPSF{ano_00_k35_psi.eps}



2525 
\end{minipage}



2526 
\\



2527 
\caption{\it Pressure, temperature, wind vectors and streamfunction (from upperleft to lowerright)



2528 
associated with the quasigeostrophic PV anomaly. The horizontal crosssection is at level 7\,km. For clarity,



2529 
some isolines of quasigeostrophic PV are also shown.}



2530 
\label{anofield}



2531 
\end{center}



2532 
\end{figure}



2533 



2534 
\vspace{0.3cm}



2535 
\begin{center}



2536 
\begin{tabular}{lll}



2537 
\hline



2538 
STREAM & Streamfunction ($m^2/s$) & ANO\_20060116\_18 \\



2539 
THETA & Potential temperature ($K$) & ANO\_20060116\_18 \\



2540 
U & Velocity component in X direction ($m/s$) & ANO\_20060116\_18 \\



2541 
V & Velocity in Y direction ($m/s$) & ANO\_20060116\_18 \\



2542 
T & Temperature ($^\circ C$) & ANO\_20060116\_18 \\



2543 
P & Pressure ($hPa$) & ANO\_20060116\_18 \\



2544 
\hline



2545 
\end{tabular}



2546 
\end{center}



2547 



2548 
\vspace{0.3cm}



2549 
\noindent



2550 
Finally, the question emerges whether the convergence of the algorithm was sufficient. In order to assess the



2551 
quality of convergence a diagnostic tool is presented in section~5.2, which allows to calculate the



2552 
quasigeostrophic PV according to its definition and then to compare it with the anomaly given at



2553 
the beginning of the PV inversion.



2554 



2555 



2556 
\subsubsection{Preparing the next step for iteration [pvin4]}



2557 



2558 



2559 
\noindent



2560 
At this point, the quasigeostrophic PV anomaly is inverted and the associated streamfunction and



2561 
other flow fields have been determined. But the inversion problem is not yet solved. In fact,



2562 
we specified an anomaly of Ertel's PV, but in the previous step a anomaly of quasigeostrophic PV



2563 
was inverted. By the above reasoning, we expect the two anomalies to be closely linked. However,



2564 
inversion of the quasigeostrophic PV equation is a linear problem, whereas the inversion of Ertel's



2565 
PV is inherently nonlinear. That is why some further iteration must be performed.\\



2566 



2567 
\noindent



2568 
In a first step, we take the basic flow fields (temperature, velocity, pressure) from the anomaly (ANO)



2569 
file and subtract it from the MOD file (Fig.\,\ref{prepiter}). This is done by call to\\[0.3cm]



2570 
\begin{center}



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



2572 
\begin{verbatim}



2573 
inversion.sh pvin3



2574 
\end{verbatim}



2575 
\end{minipage}



2576 
\end{center}



2577 



2578 
\begin{figure}[t]



2579 
\begin{center}



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



2581 
\ForceWidth{1.0\textwidth}



2582 
\BoxedEPSF{prepiter.eps}



2583 
\end{minipage}



2584 
\\



2585 
\caption{\it Preparation for the next iteration: The fields (pressure, temperature, velocity, potential



2586 
temperature) from iteration n and the anomalies of iteration n are combined to give the fields at iteration



2587 
n+1.}



2588 
\label{prepiter}



2589 
\end{center}



2590 
\end{figure}



2591 



2592 
\vspace{0.3cm}



2593 
\noindent



2594 
There is one external parameter involved in this step. It is found in the NUMERICS section of the



2595 
parameter file and is called ALPHA. This value determines which fraction of the perturbation has to



2596 
be subtracted from the MOD fields. Formally, we do the following transformation (see also Fig.\,\ref{prepiter}):



2597 



2598 
\begin{eqnarray*}



2599 
u^{n+1}_{MOD} & \leftarrow & u^n_{MOD}  \alpha \cdot u^n_{ANO}\\[0.2cm]



2600 
v^{n+1}_{MOD} & \leftarrow & v^n_{MOD}  \alpha \cdot v^n_{ANO}\\[0.2cm]



2601 
T^{n+1}_{MOD} & \leftarrow & T^n_{MOD}  \alpha \cdot T^n_{ANO}\\[0.2cm]



2602 
p^{n+1}_{MOD} & \leftarrow & p^n_{MOD}  \alpha \cdot p^n_{ANO}



2603 
\end{eqnarray*}



2604 



2605 
\vspace{0.3cm}



2606 
\noindent



2607 
A value $\alpha < 1$ slows down the convergence because only part of the inversion result is used. On the other hand, it was found that the stability of the algorithm is enhanced by values $\alpha < 1$. In



2608 
the case study a value of 0.5 was used, and good convergence reached, whereas a value of 1 resulted in unrealistic smallscale perturbations which diverge with continuing iterative steps.\\



2609 



2610 
\noindent



2611 
The variables which are adjusted in the MOD\_20060116\_18 file are given in the following table and an example



2612 
for the adjustment is shown in Fig.\,\ref{prepsub}:



2613 



2614 
\begin{center}



2615 
\begin{tabular}{lll}



2616 
\hline



2617 
THETA & Potential temperature & MOD\_20060115\_18 \\



2618 
U & Velocity component in X direction ($m/s$) & MOD\_20060115\_18 \\



2619 
V & Velocity in Y direction ($m/s$) & MOD\_20060115\_18 \\



2620 
T & Temperature ($^\circ C$) & MOD\_20060115\_18 \\



2621 
P & Pressure ($hPa$) & MOD\_20060115\_18 \\



2622 
\hline



2623 
\end{tabular}



2624 
\end{center}



2625 



2626 
\vspace{0.5cm}



2627 
\noindent



2628 
It might be valuable to study the convergence of the inversion in greater detail. To this aim, the



2629 
master script allows the call:\\[0.3cm]



2630 
\begin{center}



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



2632 
\begin{verbatim}



2633 
inversion.sh pvin5



2634 
\end{verbatim}



2635 
\end{minipage}



2636 
\end{center}



2637 



2638 
\vspace{0.3cm}



2639 
\noindent



2640 
which makes a copy of the momentarily reached MOD and ANO file. For this save option to work, the



2641 
SAVEITER flag in the NUMERICS section of the parameter file must be set to "yes".



2642 
\begin{center}



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



2644 
\begin{verbatim}



2645 
BEGIN NUMERICS



2646 
SAVEITER = yes;



2647 
END NUMERICS



2648 
\end{verbatim}



2649 
\end{minipage}



2650 
\end{center}



2651 



2652 
\begin{figure}[t]



2653 
\begin{center}



2654 
\begin{minipage}[t]{8.3cm}



2655 
\ForceWidth{1.0\textwidth}



2656 
\BoxedEPSF{prep_k35_org_ano.eps}



2657 
\end{minipage}



2658 
\hspace{1cm}



2659 
\begin{minipage}[t]{8.3cm}



2660 
\ForceWidth{1.0\textwidth}



2661 
\BoxedEPSF{prep_k35_mod.eps}



2662 
\end{minipage}



2663 
\\



2664 
\caption{\it Left: Temperature in the modified MOD file (in color) and in the anomaly ANO file (with contour lines)



2665 
at the end of the quasigeostrophic PV inversion. The temperature anomaly in ANO is subtracted from the initial temperature in MOD. Right: New modified temperature in the MOD file after the subtraction.



2666 
This new modified field gives the next entry for the PV inversion. }



2667 
\label{prepsub}



2668 
\end{center}



2669 
\end{figure}



2670 



2671 



2672 
\subsubsection{Convergence after several iterations [pvin5]}



2673 



2674 
At this stage, the next iterative steps can be performed. The MOD file has been adjusted according to



2675 
the outcome of the quasigeostrophic PV inversion. The modified wind and temperature fields can be used to calculate a new ErtelPV field, which in turn can be compared with the prespecified aimPV field.



2676 
Hopefully, after several iterative steps, the recurrently calculated PV of the MOD file converges



2677 
toward the aimPV. If so, the nonlinear inversion problem for Ertel's PV is solved. Figure\,\ref{convergence}



2678 
illustrates the convergence for the case study. Note that the iterative steps need not be



2679 
launched by hand, but can be run automatically with call to\\[0.3cm]



2680 
\begin{center}



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



2682 
\begin{verbatim}



2683 
inversion.sh pvin



2684 
\end{verbatim}



2685 
\end{minipage}



2686 
\end{center}



2687 



2688 
\vspace{0.3cm}



2689 
\noindent



2690 
The script will run through all inversion steps and will loop through the number of iterations specified in NOFITER in the NUMERICS section of the parameter file.



2691 
\begin{center}



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



2693 
\begin{verbatim}



2694 
BEGIN NUMERICS



2695 
NOFITER = 6;



2696 
END NUMERICS



2697 
\end{verbatim}



2698 
\end{minipage}



2699 
\end{center}



2700 



2701 
\vspace{0.3cm}



2702 
\noindent



2703 
The panels of Fig.\,\ref{convergence} show the streamfunction at three consecutive iterative steps.



2704 
Note that the amplitiude of the streamfunction clearly decreases, indicating that the iteration



2705 
process converges.



2706 
Generally, after six steps good convergence is reached. This is further supported by



2707 
the standard output during the quasigeostrophic inversion (see listing on page~38). The final line of this output



2708 
is shown below for several consecutive iterations (the first one corresponding to the line on page~38):



2709 



2710 
\begin{figure}[t]



2711 
\begin{center}



2712 
\begin{minipage}[t]{5cm}



2713 
\ForceWidth{1.0\textwidth}



2714 
\BoxedEPSF{con_01_k35_psi.eps}



2715 
\end{minipage}



2716 
\hspace{0cm}



2717 
\begin{minipage}[t]{5cm}



2718 
\ForceWidth{1.0\textwidth}



2719 
\BoxedEPSF{con_02_k35_psi.eps}



2720 
\end{minipage}



2721 
\hspace{0cm}



2722 
\begin{minipage}[t]{5cm}



2723 
\ForceWidth{1.0\textwidth}



2724 
\BoxedEPSF{con_03_k35_psi.eps}



2725 
\end{minipage}



2726 
\\



2727 
\caption{\it streamfunction for three iterative steps. Note the different contour intervals



2728 
for the three sub panels. The decrease of the amplitude clearly indicates a convergence of the



2729 
PV inversion.}



2730 
\label{convergence}



2731 
\end{center}



2732 
\end{figure}



2733 



2734 
\begin{center}



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



2736 
\begin{verbatim}



2737 
psigauge = 0.309E+03 deltasq = 0.180E10



2738 
psigauge = 0.128E+03 deltasq = 0.548E11



2739 
psigauge = 0.570E+02 deltasq = 0.233E11



2740 
psigauge = 0.265E+02 deltasq = 0.117E11



2741 
\end{verbatim}



2742 
\end{minipage}



2743 
\end{center}



2744 



2745 
\vspace{0.3cm}



2746 
\noindent



2747 
Note that with each iterative step "less" streamfunction is produced. Since the basic perturbation



2748 
fields (temperature, pressure, horizontal wind, potential temperature) are directly proportional



2749 
to the streamfunction, the corrections to the flow fields also



2750 
become smaller and smaller.



2751 



2752 
\subsection{Post processing  Changing the P file [post]}



2753 



2754 
The PV inversion was completed in the previous step and the modified temperature, pressure and wind field is available in the run files. It is the aim of the post processing steps to bring these modified



2755 
fields back to the format of the input fields. The steps are summarised in Fig.\,\ref{postproc}.\\



2756 



2757 
\begin{figure}[t]



2758 
\begin{center}



2759 
\begin{minipage}[t]{12cm}



2760 
\ForceWidth{1.0\textwidth}



2761 
\BoxedEPSF{postproc.eps}



2762 
\end{minipage}



2763 
\\



2764 
\caption{\it Different steps of post processing. Additionally, to the left of the different



2765 
blocks the main flow of information is indicated. For instance, MOD $\rightarrow$ GEO means



2766 
that the input is taken from the MOD file and the output written to the GEO file. }



2767 
\label{postproc}



2768 
\end{center}



2769 
\end{figure}



2770 



2771 
\noindent



2772 
But at first the input field and the link to the modified output field must be made available with a call to:\\[0.3cm]



2773 
\begin{center}



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



2775 
\begin{verbatim}



2776 
inversion.sh post1



2777 
\end{verbatim}



2778 
\end{minipage}



2779 
\end{center}



2780 



2781 
\vspace{0.3cm}



2782 
\noindent



2783 
The P file is copied from the input directory to the output directory and a symbolic link is set to the



2784 
MOD file in the run directory.



2785 



2786 



2787 



2788 



2789 
\subsubsection{Rotating to the geographical latitude/longitude grid [post2]}



2790 



2791 
Figure~\ref{backrot} illustrates the step which is necessary to transform the meteorological fields to a geographical coordinate system. The left panel shows the temperature at 2000\,m above sea level in the localcartesian coordinate system, the right panel is the same field, but projected back to a geographical grid. Note that in this latter frame only a small section is filled by temperature values, whereas most are



2792 
undefined.\\



2793 



2794 
\begin{figure}[t]



2795 
\begin{center}



2796 
\begin{minipage}[t]{8.3cm}



2797 
\ForceWidth{1.0\textwidth}



2798 
\BoxedEPSF{rot_k35_rot_t.eps}



2799 
\end{minipage}



2800 
\hspace{1.2cm}



2801 
\begin{minipage}[t]{8.3cm}



2802 
\ForceWidth{1.0\textwidth}



2803 
\BoxedEPSF{rot_k35_geo_t.eps}



2804 
\end{minipage}



2805 
\\



2806 
\caption{\it Modified temperature field at 7\,km above ground in the rotated coordinate system (left) and transformed back into the geographical latitude/longitude grid (right). Additionally



2807 
the wind vectors are shown in both coordinate systems. Note that the transformation of the vectorial



2808 
wind is much more complicated than the transformation of the scalar temperature.}



2809 
\label{backrot}



2810 
\end{center}



2811 
\end{figure}



2812 



2813 
\noindent



2814 
The rotation itself is the reverse transformation performed in the preparatory steps (section~4.4.2). It is called by (and its effect shown in Fig.\,\ref{backrot})\\[0.3cm]



2815 
\begin{center}



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



2817 
\begin{verbatim}



2818 
inversion.sh post2



2819 
\end{verbatim}



2820 
\end{minipage}



2821 
\end{center}



2822 



2823 
\vspace{0.3cm}



2824 
\noindent



2825 
and its behaviour is controlled by the following entries in the parameter file:



2826 
\begin{center}



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



2828 
\begin{verbatim}



2829 
BEGIN GRID



2830 
GEO_XMIN = 180.;



2831 
GEO_NX = 361 ;



2832 
GEO_DX = 1.;



2833 
GEO_YMIN = 0.;



2834 
GEO_NY = 91 ;



2835 
GEO_DY = 1.;



2836 
CLON = 65.;



2837 
CLAT = 45.;



2838 
CROT = 0.;



2839 
END GRID



2840 
\end{verbatim}



2841 
\end{minipage}



2842 
\end{center}



2843 



2844 
\vspace{0.3cm}



2845 
\noindent



2846 
The first six parameters define the geographical grid of the input files. More specifically the left grid boundary is at $180^\circ$\,W (GEO\_XMIN), the number of grid points in the zonal (westtoeast) direction is 361 (GEO\_NX) and



2847 
the grid resolution in zonal direction is $1^\circ$\,longitude (GEO\_DX). Analogously, the parameters of the



2848 
grid in meridional (southtonorth) direction are given: Lowest latitude $0^\circ$\,N (GEO\_YMIN), number of grid



2849 
points 91 (GEO\_NY) and grid resolution $1^\circ$\,latitude (GEO\_DY). The rotation itself is described by



2850 
the last three entries (CLON, CLAT and CROT), which have the exactly same meaning as given in section~4.4.2\\



2851 



2852 
\noindent



2853 
Note again, that the transformation of a scalar quantity is considerably easier than the transformation of a vectorial quantity. In the latter case, the specific transformation



2854 
of the vectorial components must be correctly handled. The transformation algorithm for a scalar quantity is based



2855 
upon the transformation rules given in Appendix~9.1 (Fortran subroutines {\em lmtolms} and {\em phtophs}). Essentially, the



2856 
steps are the reverse of the ones presented in section~4.4.2.



2857 
The rotated (quasicartesian) latitude/longitude corresponding to a grid point in the geographical coordinate system is obtained in the following way. Firstly, the geographical latitude/longitude $\phi_{geo},\lambda_{geo}$ is transformed into



2858 



2859 
\begin{eqnarray*}



2860 
\lambda'_{rot} & = & lmtolms(\phi_{geo},\lambda_{geo},90\phi_{cen},\lambda_{cen}180) \\[0.2cm]



2861 
\phi'_{rot} & = & phtophs(\phi_{geo},\lambda_{geo},90\phi_{cen},\lambda_{cen}180))



2862 
\end{eqnarray*}



2863 



2864 
\vspace{0.3cm}



2865 
\noindent



2866 
where the centre of the PV anomaly (the centre of the quasicartesian coordinate system in geographical



2867 
latitude/longitude coordinates) is given by $\phi_{cen}$ and $\lambda_{cen}$ (CLAT and CLON in the



2868 
parameter file, see above).



2869 
These two



2870 
values are then, in a second rotation, transformed into rotated latitude and longitude:



2871 



2872 
\begin{eqnarray*}



2873 
\lambda_{rot} & = & 90+lmtolms(\phi'_{rot},\lambda'_{rot}90,90+\alpha,180) \\[0.2cm]



2874 
\phi_{rot} & = & phtophs(\phi'_{rot},\lambda'_{rot}90,90+\alpha,180)



2875 
\end{eqnarray*}



2876 



2877 
\vspace{0.3cm}



2878 
\noindent



2879 
where $\alpha$ is the rotation angle given in the parameter file (CROT, see below). With the correspondence



2880 
$(\lambda_{geo},\phi_{geo}) \leftrightarrow (\lambda_{rot},\phi_{rot})$ it is straightforward to perform



2881 
the transformation of a scalar field. In the program package this is done by means of linear interpolation.\\



2882 



2883 
\noindent



2884 
The fields which are written to GEO\_20060116\_18 are given in the following table:\\



2885 



2886 
\begin{center}



2887 
\begin{tabular}{lll}



2888 
\hline



2889 
U & Zonal velocity ($m/s$) & GEO\_20060116\_18 \\



2890 
V & Meridional velocity ($m/s$) & GEO\_20060116\_18 \\



2891 
T & Temperature ($^\circ C$) & GEO\_20060116\_18 \\



2892 
P & Pressure ($hPa$) & GEO\_20060116\_18 \\



2893 
\hline



2894 
\end{tabular}



2895 
\end{center}



2896 
\vspace{0.5cm}



2897 



2898 
\begin{figure}[t]



2899 
\vspace{3cm}



2900 
\begin{center}



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



2902 
\ForceWidth{1.0\textwidth}



2903 
\BoxedEPSF{weight1.eps}



2904 
\end{minipage}



2905 
\\[4.2cm]



2906 
\caption{\it Weighting function (in color) for the insertion of the modified fields into the original



2907 
P file. Additionally, the distance from the boundary of insertion region is plotted as contour lines (positive



2908 
inside the region, negative outside region). The zero contour line of the distance corresponds to the



2909 
boundary of the insertion domain.}



2910 
\label{comparep}



2911 
\end{center}



2912 
\end{figure}



2913 



2914 



2915 
\subsubsection{Transformation from height to hybrid ECMWF coordinates [post3]}



2916 



2917 
\begin{figure}[t]



2918 
\vspace{1cm}



2919 
\begin{center}



2920 
\begin{minipage}[t]{8.2cm}



2921 
\ForceWidth{1.0\textwidth}



2922 
\BoxedEPSF{inp_tt_k38.eps}



2923 
\end{minipage}



2924 
\hspace{0.7cm}



2925 
\begin{minipage}[t]{8.2cm}



2926 
\ForceWidth{1.0\textwidth}



2927 
\BoxedEPSF{out_tt_k38.eps}



2928 
\end{minipage}



2929 
\\[2.8cm]



2930 
\begin{minipage}[t]{8.2cm}



2931 
\ForceWidth{1.0\textwidth}



2932 
\BoxedEPSF{inp_uv_k38.eps}



2933 
\end{minipage}



2934 
\hspace{0.7cm}



2935 
\begin{minipage}[t]{8.2cm}



2936 
\ForceWidth{1.0\textwidth}



2937 
\BoxedEPSF{out_uv_k38.eps}



2938 
\end{minipage}



2939 
\\[1.5cm]



2940 
\caption{\it Comparison of original (with PV streamer, left) and modified (without PV streamer, right)



2941 
temperature (top panels) and velocity/wind vectors bottom panels) at model level 38 of the ECMWF grid.}



2942 
\label{comparep2}



2943 
\end{center}



2944 
\end{figure}



2945 



2946 
\noindent



2947 
The vertical coordinate in the backrotated field is still geometrical height, whereas the input



2948 
fields are given on the hybrid ECMWF grid. In this step the transformation to the hybrid vertical coordinate is performed. The call is:\\[0.3cm]



2949 
\begin{center}



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



2951 
\begin{verbatim}



2952 
inversion.sh post3



2953 
\end{verbatim}



2954 
\end{minipage}



2955 
\end{center}



2956 



2957 
\vspace{0.3cm}



2958 
\noindent



2959 
The transformation is based upon a cubic spline interpolation. The fields which are overwritten



2960 
in the P20060116\_18 file are:\\



2961 



2962 
\begin{center}



2963 
\begin{tabular}{lll}



2964 
\hline



2965 
U & Velocity component in zonal direction ($m/s$) & P20060115\_18 \\



2966 
V & Velocity in meridional direction ($m/s$) & P20060115\_18 \\



2967 
T & Temperature ($^\circ C$) & P20060115\_18 \\



2968 
PS & Surface pressure ($hPa$) & P20060115\_18 \\



2969 
\hline



2970 
\end{tabular}



2971 
\end{center}



2972 



2973 
\vspace{0.3cm}



2974 
\noindent



2975 
An impression of the modified flow field is given in Fig.\,\ref{comparep2}.



2976 



2977 



2978 
\subsubsection{Calculating secondary fields [post4]}



2979 



2980 
\begin{figure}[t]



2981 
\begin{center}



2982 
\begin{minipage}[t]{7.8cm}



2983 
\ForceWidth{1.0\textwidth}



2984 
\BoxedEPSF{inp_pv.th_45n.eps}



2985 
\end{minipage}



2986 
\hspace{0cm}



2987 
\begin{minipage}[t]{7.8cm}



2988 
\ForceWidth{1.0\textwidth}



2989 
\BoxedEPSF{out_pv.th_45n.eps}



2990 
\end{minipage}



2991 
\\[1.5cm]



2992 
\begin{minipage}[t]{7.8cm}



2993 
\ForceWidth{1.0\textwidth}



2994 
\BoxedEPSF{inp_pv_350hPa.eps}



2995 
\end{minipage}



2996 
\hspace{0cm}



2997 
\begin{minipage}[t]{7.8cm}



2998 
\ForceWidth{1.0\textwidth}



2999 
\BoxedEPSF{out_pv_350hPa.eps}



3000 
\end{minipage}



3001 
\\[1.5cm]



3002 
\caption{\it Comparison of original (with PV streamer, left) and modified (without PV streamer, right)



3003 
potential temperature and PV in a west/east cross section along the $45^\circ$\,N latitude circle (top panels)



3004 
and PV at 350\,hPa (bottom panels)}



3005 
\label{compares}



3006 
\end{center}



3007 
\end{figure}



3008 



3009 
In this last step, Ertel's PV and potential temperature are calculated on the ECMWF grid. This allows



3010 
to directly check whether the PV anomaly was removed and to investigate how the corresponding potential temperature is



3011 
changed in the absence of the PV anomaly. The call to this step is:\\[0.3cm]



3012 
\begin{center}



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



3014 
\begin{verbatim}



3015 
inversion.sh post4



3016 
\end{verbatim}



3017 
\end{minipage}



3018 
\end{center}



3019 



3020 
\vspace{0.3cm}



3021 
\noindent



3022 
The program creates the S file and writes the following two fields onto it:\\[0.3cm]



3023 



3024 
\begin{center}



3025 
\begin{tabular}{lll}



3026 
\hline



3027 
TH & Potential temperature ($K$) & S20060116\_18 \\



3028 
PV & Ertel's PV ($pvu$) & S20060116\_18 \\



3029 
\hline



3030 
\end{tabular}



3031 
\end{center}



3032 



3033 
\vspace{0.3cm}



3034 
\noindent



3035 
Figure~\ref{compares} shows a direct comparison of the input and of the output ErtelPV and potential temperature.



3036 
It can clearly be seen that the stratospheric PV streamer is essentially removed from the output



3037 
file. Moreover, the potential temperature is considerably changed. Whereas in the presence of the



3038 
PV streamer the isolines of potential temperature are pulled upwards, they are much more horizontally



3039 
aligned in its absence. This "vacuum cleaner" effect of upperlevel PV structures is well documented



3040 
and predicted by theoretical considerations (see also Fig.\,2 in section~1).



3041 



3042 



3043 
% 



3044 
% PV Inversion for Idealised Cases



3045 
% 



3046 



3047 
\newpage



3048 



3049 



3050 
\section{Diagnostic Tools}



3051 



3052 
\subsection{A consistency check for the boundary condition [diag1]}



3053 



3054 
The problem posed by the inversion equations constitutes a von Neumann boundary value problem.



3055 
It has only a solution if some additional conditions are fulfilled between the interior PV distribution and the values of potential temperature on the upper and lower boundaries and the values of horizontal wind components on the lateral sides of the inversion domain. What follows is a



3056 
detailed description of these additional requirements. The discussion very closely follows the one given



3057 
in Fehlmann (1997).\\



3058 



3059 
\noindent



3060 
To facilitate the derivation, let's define a new vertical coordinate such theta the quasigeostrophic PV q can be expressed as the divergence of a vector field $\vec{E}$. Let\\[0.2cm]



3061 



3062 
\[



3063 
\eta(z) =  \frac{p_0(z)}{g f^2} \quad \quad \mbox{and} \quad \quad



3064 
\Delta_h = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}



3065 
\]



3066 



3067 
\vspace{0.3cm}



3068 
\noindent



3069 
where $p_o(z)$ is the reference profile of pressure, $f$ is the Coriolis parameter and $g$ is the Earth's gravity. Note that $\eta$ is a scaled pressure coordinate. If we use the hydrostatic assumption $\partial{\eta}/{\partial z} = \rho_0/f^2$, the following expression for the quasigeostrophic PV q results:



3070 



3071 
\[



3072 
q = \Delta_h\psi + \frac{\partial}{\partial \eta} ( \frac{\rho_0^2}{f^2 N_0^2}



3073 
\cdot \frac{\partial \psi}{\partial \eta} )



3074 
\]



3075 



3076 
\vspace{0.3cm}



3077 
\noindent



3078 
Here, $\rho_o(z)$ and $N_o^2(z)$ are the reference profiles of density and squared BruntVais\"al\"a frequency,



3079 
respectively, and $\psi$ is the streamfunction. Hence, q is now expressed as the divergence of a vector field ${\vec E}$, and Gauss' theorem can be



3080 
applied:



3081 



3082 
\[



3083 
\int_G q(x,y,\eta) dx dy d\eta = \int_{\partial G} \vec{E}(x,y,\eta) \vec{d\sigma} \quad \quad



3084 
\mbox{with} \quad \quad



3085 
\vec{E} = ( \frac{\partial \psi}{\partial x}, \frac{\partial \psi}{\partial y},



3086 
\frac{\rho_0^2}{f^2 N_0^2} \cdot \frac{\partial \psi}{\partial \eta} )



3087 
\]



3088 



3089 
\vspace{0.3cm}



3090 
\noindent



3091 
Note that this integration has to be carried out in the coordinate system (x,y,$\eta$). If we now specify



3092 
the inversion domain B as a twodimensional domain in the (x,y)plane and $[z_a,z_b]$ be the vertical domain, then $G=B \times [z_a,z_b]$ in the Cartesian space. This compatibility condition can be



3093 
transformed back into the Cartesian coordinate system and the integration be split over the boundary of



3094 
the domain G into an integration over the surface, the lid and the lateral boundaries of the domain. It then



3095 
follows that



3096 



3097 
\begin{eqnarray*}



3098 
\int_G \rho_0(z) q(x,y,z) dx dy dz & = &



3099 
\int _{z_a}^{z_b} (\int_{\partial B} \vec{v}(x,y,z) \cdot \vec{dr} ) \rho_0(z) dz\\[4mm]



3100 
&  & \int_B \frac{\rho_0(z_a) g f \theta^*(x,y,z_a)}{\theta_0(z_a) N_0^2(z_a)} dx dy \\[4mm]



3101 
& + & \int_B \frac{\rho_0(z_b) g f \theta^*(x,y,z_b)}{\theta_0(z_b) N_0^2(z_b)} dx dy \\



3102 
\end{eqnarray*}



3103 



3104 
\vspace{0.3cm}



3105 
\noindent



3106 
Only if this condition is fulfilled, does there exist a unique solution (up to a constant). This means



3107 
that not every formulation of the inversion problem needs to have a solution.\\



3108 



3109 
\noindent



3110 
The inversion program has to check whether the just discussed integral constraint is fulfilled to



3111 
satisfactory accuracy. In order to estimate, how big these inconsistencies really are, we ask for a constant perturbation potential temperature $\theta^*$ which is needed in order to fulfill the integral



3112 
constraint. If readily follows from the above discussion that this value is:



3113 



3114 
\begin{eqnarray*}



3115 
\theta^* & = &



3116 
( \int_G \rho_0(z) q(x,y,z) dx dy dz 



3117 
\int _{z_a}^{z_b} (\int_{\partial B} \vec{v}(x,y,z) \cdot \vec{dr} ) \rho_0(z) dz )\\[4mm]



3118 
& &



3119 
(



3120 
 \int_B \frac{\rho_0(z_a) g f}{\theta_0(z_a) N_0^2(z_a)} dx dy



3121 
+ \int_B \frac{\rho_0(z_b) g f}{\theta_0(z_b) N_0^2(z_b)} dx dy )^{1}\\



3122 
\end{eqnarray*}



3123 



3124 
\vspace{0.3cm}



3125 
\noindent



3126 
If this quantity is calculated for the present example, $\theta^*=$2.8\,K results. This value seems



3127 
sufficiently small that problem can be solved despite the inconsistency. But note that the problems due



3128 
to the illposedness of the problem are not well understood (see references in Fehlmann, 1997).\\



3129 



3130 
\noindent



3131 
The program handling the consistency check related to the boundary values is:\\[0.3cm]



3132 
\begin{center}



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



3134 
\begin{verbatim}



3135 
inversion.sh diag1



3136 
\end{verbatim}



3137 
\end{minipage}



3138 
\end{center}



3139 



3140 
\vspace{0.3cm}



3141 
\noindent



3142 
The output from this program is given in the following table (the labels A,B,C and D are entered for easier



3143 
reference in the text):\\[0.3cm]



3144 
\begin{center}



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



3146 
\begin{verbatim}



3147 
A integ = 1.7558686E+12



3148 
B denombot = 5.6644076E+11



3149 
C denomtop = 5.9614413E+09



3150 
D denom = 5.6047934E+11



3151 
theta adjustment = 3.132798



3152 
theta shift @ top = 0.000000 613.2715



3153 
theta shift @ bot = 0.000000 271.9355



3154 
\end{verbatim}



3155 
\end{minipage}



3156 
\end{center}



3157 



3158 
\vspace{0.3cm}



3159 
\noindent



3160 
The first line is the integral which includes the quasigeostrophic PV in the interior, the



3161 
velocity integrals on the laterals boundaries and the potential temperature integrals at the



3162 
lower and upper boundary:



3163 



3164 
\begin{eqnarray*}



3165 
A & = & \int_G \rho_0(z) q(x,y,z) dx dy dz 



3166 
\int _{z_a}^{z_b} (\int_{\partial B} \vec{v}(x,y,z) \cdot \vec{dr} ) \rho_0(z) dz\\[6mm]



3167 
& & + \int_B \frac{\rho_0(z_a) g f \theta^*(x,y,z_a)}{\theta_0(z_a) N_0^2(z_a)} dx dy



3168 
 \int_B \frac{\rho_0(z_b) g f \theta^*(x,y,z_b)}{\theta_0(z_b) N_0^2(z_b)} dx dy \\



3169 
\end{eqnarray*}



3170 



3171 
\vspace{0.3cm}



3172 
\noindent



3173 
In exact fulfillment of the integral constraint, this term A would vanish, i.e. $A=0$.



3174 
The second and third line give the integrals over the upper and lower boundary only:



3175 



3176 
\begin{eqnarray*}



3177 
B & = & \int_B \frac{\rho_0(z_a) g f \theta^*(x,y,z_a)}{\theta_0(z_a) N_0^2(z_a)} dx dy\\[6mm]



3178 
C & = & \int_B \frac{\rho_0(z_b) g f \theta^*(x,y,z_b)}{\theta_0(z_b) N_0^2(z_b)} dx dy\\



3179 
\end{eqnarray*}



3180 



3181 
\vspace{0.3cm}



3182 
\noindent



3183 
The difference between these two lines is the contents of the next line:



3184 



3185 
\begin{eqnarray*}



3186 
D & = & \int_B \frac{\rho_0(z_a) g f \theta^*(x,y,z_a)}{\theta_0(z_a) N_0^2(z_a)} dx dy



3187 
 \int_B \frac{\rho_0(z_b) g f \theta^*(x,y,z_b)}{\theta_0(z_b) N_0^2(z_b)} dx dy\\



3188 
\end{eqnarray*}



3189 



3190 
\noindent



3191 
Note that this term D corresponds to the last two terms in expression~A.



3192 
With these expressions an estimate can be made for the correction which would eliminate the inconsistency



3193 
of the boundary conditions with respect to the interior distribution of quasigeostrophic PV. This adjustment



3194 
of potential temperature is the next line



3195 



3196 
\begin{eqnarray*}



3197 
\theta^* & = & A / D



3198 
\end{eqnarray*}



3199 



3200 
\vspace{0.3cm}



3201 
\noindent



3202 
and corresponds with the aforementioned expression for $\theta^*$.



3203 
Finally, in the last two lines the shift in potential temperature at the upper and at the



3204 
lower boundary is given, together with the mean potential temperatures at these two boundaries.



3205 
At the moment, no shift (value 0 in the output lines) is performed, i.e. the inconsistency in the boundary conditions



3206 
is not "solved".



3207 



3208 



3209 
\subsection{Convergence of the quasigeostrophic inversion [diag2]}



3210 



3211 
The key numerical algorithm is the inversion of the quasigeostrophic PV equation (see section~4.5.3).



3212 
It relates the PV to the streamfunction by\\[0.3cm]



3213 



3214 
\[



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



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



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



3218 
\]



3219 



3220 
\vspace{0.3cm}



3221 
\noindent



3222 
with the boundary conditions of potential temperature at the lower and upper lid, and of horizontal wind



3223 
components at the lateral sides of the inversion domain:\\[0.3cm]



3224 



3225 
\[



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



3227 
\quad \quad



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



3229 
\quad \quad



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



3231 
\]



3232 



3233 
\begin{figure}[t]



3234 
\begin{center}



3235 
\begin{minipage}[t]{8.3cm}



3236 
\ForceWidth{1.0\textwidth}



3237 
\BoxedEPSF{qgpv_k35_calc.eps}



3238 
\end{minipage}



3239 
\hspace{1cm}



3240 
\begin{minipage}[t]{8.3cm}



3241 
\ForceWidth{1.0\textwidth}



3242 
\BoxedEPSF{qgpv_k35_spec.eps}



3243 
\end{minipage}



3244 
\\



3245 
\caption{\it Comparison of calculated (left) and predescribed (right) quasigeostrophic PV. }



3246 
\label{compqg}



3247 
\end{center}



3248 
\end{figure}



3249 



3250 
\vspace{0.3cm}



3251 
\noindent



3252 
The numerical solution by means of the successive overrelaxation (SOR) technique



3253 
then completely determines the perturbations of potential temperature and of horizontal wind. These meteorological fields can be expressed



3254 
in the following way by the streamfunction (see also section~2):



3255 



3256 
\begin{eqnarray*}



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



3258 
u & = &  \frac{\partial \psi}{\partial y}\\



3259 
v & = & \frac{\partial \psi}{\partial x}\\



3260 
\end{eqnarray*}



3261 



3262 
\vspace{0.3cm}



3263 
\noindent



3264 
This allows a simple test whether the inversion of the quasigeostrophic equation was successful. Indeed, the



3265 
quasigeostrophic PV can be calculated from the perturbation fields:



3266 



3267 
\[



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



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



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



3271 
\]



3272 



3273 
\vspace{0.3cm}



3274 
\noindent



3275 
If the inversion was successful, this calculated quasigeostrophic PV ($q_c$) must coincide with the



3276 
predescribed PV ($q_a$). For the present case study, the validation is shown in Fig.\,\ref{compqg} for the first step



3277 
of the iteration. The two panels show the calculated and the predescribed PV, respectively, and it becomes



3278 
immediately clear that a reasonable degree of convergence was reached.\\



3279 



3280 



3281 
\noindent



3282 
The calculation of the quasigeostrophic PV is accomplished with the call\\[0.3cm]



3283 
\begin{center}



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



3285 
\begin{verbatim}



3286 
inversion.sh diag2 ANO



3287 
\end{verbatim}



3288 
\end{minipage}



3289 
\end{center}



3290 



3291 
\vspace{0.3cm}



3292 
\noindent



3293 
where the second argument (ANO) specifies that the ANO\_20060116\_18 file is taken for the



3294 
calculation. It is implictely assumed that the reference file is available on the corresponding REF file. The output of the calculation is then written to the same ANO\_20060116\_18 file,



3295 
with the name QGPV\_DIAG. Hence the table of modified fields and files is:



3296 



3297 
\begin{center}



3298 
\begin{tabular}{lll}



3299 
\hline



3300 
QGPV\_DIAG & Calculated quasigeostrophic PV ($1/s$) & ANO\_20060116\_18 \\



3301 
\hline



3302 
\end{tabular}



3303 
\end{center}



3304 
\vspace{0.3cm}



3305 



3306 
\subsection{Difference between two files [diag3]}



3307 



3308 
\begin{figure}[t]



3309 
\begin{center}



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



3311 
\vspace{1cm}



3312 
\ForceWidth{1.0\textwidth}



3313 
\BoxedEPSF{difference.eps}



3314 
\end{minipage}



3315 
\\[0.5cm]



3316 
\caption{\it Difference between the ORG and the MOD file in a vertical cross section. In color, the difference of the meridional



3317 
wind velocity is shown. The difference of potential temperature is given by the thin lines (black:positive,



3318 
blue:negative). The bold black line is the 2.5 isoline of the difference of Ertel's PV. }



3319 
\label{diff}



3320 
\end{center}



3321 
\end{figure}



3322 



3323 
The impact of a PV anomaly is most easily determined if the meteorological fields (temperature, velocity,



3324 
pressure,...) in the presence of the anomaly is subtracted from the corresponding fields in the absence



3325 
of the anomaly. The change in the meteorological fields can then be attributed to the PV anomaly. In order



3326 
to facilitate this with/without comparison, a diagnostic tool is provide which calculates the difference.



3327 
For instance, the call to



3328 
\begin{center}



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



3330 
\begin{verbatim}



3331 
inversion.sh diag3 ORG MOD



3332 
\end{verbatim}



3333 
\end{minipage}



3334 
\end{center}



3335 



3336 
\vspace{0.3cm}



3337 
\noindent



3338 
takes the ORG and MOD file in the run directory and calculates the difference of all fields which



3339 
are available on both input files. The difference fields are then written to a new file with prefix DIA.



3340 
An example is shown in Fig.\,\ref{diff}.



3341 



3342 
% 



3343 
% Final Remarks



3344 
% 



3345 



3346 
\newpage



3347 



3348 
\section{Final Remarks and Outlook}



3349 



3350 
Several additional extensions are possible for the PV inversion tool introduced in this study. Most of them will attract a lot of interest in the research and education community. Some extensions which are and will be undertaken recently in our institute are:



3351 



3352 
\begin{itemize}



3353 
\item [a] The PV inversion is not only of interest in realcase studies, as was described in this work. Considerable insight can also be gained by looking at highly idealised experiments. Such experiments enforce the existence of a tool which allows to set set up them. In the course of this diploma thesis, such a tool was developed. However, it was considered to be preferable to limit



3354 
the thesis's documentation itself to the realcase experiments.\\



3355 
\item [b] Removing and adding, or more generally modifying the initial PV distribution, is particularly interesting if the modified state of the atmosphere is used as the input initial



3356 
and boundary values for a regional weather prediction model. This kind of experiments allows to



3357 
assess in a quantitative way how the PV field influences the future weather development. In parallel to this thesis a program package has been developed at the Institute for Atmospheric and Climate Science at ETH Zurich (IACETH) which allows to transform



3358 
the (modified) P files of the PV inversion directly into the input and boundary files for the



3359 
Climate HRM model. It is the aim that some case studies will be undertaken in this way by a



3360 
diploma student at our institute.\\



3361 
\item [c] From the quasigeostrophic omega equation it is well known that PV structures are able to



3362 
enforce vertical motions. These are particularly interesting because vertical motion leads to



3363 
condensation or convective instability, if the vertical temperature and humidity profile are suitable. Therefore,



3364 
it would be of great interest to see how the modified PV distribution is associated with different



3365 
vertical motions. Technically this means that an additional equation, the aforementioned



3366 
quasigeostrophic omega equation, has to be solved. In the coming summer semester, a student with focus on scientific



3367 
programming will attack this problem.\\



3368 
\item [d] Teaching atmospheric dynamics at the graduate level involves several challenges. One



3369 
particular challenge is associated with the abstract concept of potential vorticity. An easytouse



3370 
PV inversion tool might be of great relieve to the teacher since he can easily illustrate and



3371 
apply the concepts to real and idealised cases. It is quite certain that such a visualisation help



3372 
would be much appreciated by graduate students in atmospheric dynamics. It is one aim of the



3373 
dynamic's group at {\it IAC}ETHto improve teaching in this respect, while keeping the mathematical and physical



3374 
level of the taught contents.



3375 
\end{itemize}



3376 



3377 



3378 
\newpage



3379 
\section{Acknowledgment}



3380 



3381 



3382 
There are many persons which need to be mentioned in this acknowledgment. Perhaps, most crucially is



3383 
a heartily thanks to Rene Fehlmann who initially developed the inversion tool. I had the joy to start my atmospheric research under Rene's supervision, and his enthusiasm for mathematics, numerics and atmosphere



3384 
always was and still is very inspiring for me. \\



3385 



3386 
\noindent



3387 
Recently I read that success in science cannot be attributed solely to one's own intelligence an



3388 
"superiority". And with my long experience in science, I fully agree with that statement. Success needs



3389 
luck, and last but not least it needs excellent teachers, supervisors and colleagues. I had the great pleasure



3390 
to work together and get help from excellent scientists and teachers: Christoph Sch\"ar, Heini Wernli and Huw Davies. I appreciated very much their scientific qualities, but probably still more their warmth and



3391 
friendliness as colleagues.\\



3392 



3393 
\noindent



3394 
Many research colleagues were indirectly important in the development of this PV inversion tool. Conny Schwierz, Olivia Martius, Mischa CrociMaspoli, Sarah Kew. Certainly, thanks should also go to Linda Schlemmer who



3395 
accepted a diploma thesis based on this tool, and is now testing all aspects of it.\\



3396 



3397 
\noindent



3398 
Finally, this work was done in the frame of the postgraduate course in Computer Science at the Fernfachhochschule Schweiz. Thanks goes to all teachers, and in particular to Lars Kruse who accepted to supervise this diploma thesis and discussed with me the different steps. Although we had some



3399 
different views what should be included into this thesis ("Why do you want to send this part to the appendix? This is a diploma thesis,...") I enjoyed working with him, particularly because of his interest in the atmospheric sciences.



3400 



3401 



3402 
\newpage



3403 
\section{References}



3404 



3405 
\noindent



3406 
{\bf [1]} Acheson, D.\,J., 1990:



3407 
Elementary Fluid Dynamics.



3408 
{\sl Oxford Applied Mathematics and Computing Science Series,}, Oxford UP, 395pp.\\



3409 



3410 
\noindent



3411 
{\bf [2]} Appenzeller, C., H.~C.~Davies, and W.~A.~Norton, 1996:



3412 
Fragmentation of stratospheric intrusions.



3413 
{\sl J.~Geophys.~Res.,} {\bf 101,} 14351456.\\



3414 



3415 
\noindent



3416 
{\bf [3]} Hoskins, B.\,J., M.\,E.\,McIntyre, and A.\,W.\,Robertson, 1985:



3417 
On the use and significance of isentropic potential vorticity maps. Quart.\,J.\,Roy.\,Meteor.\,Soc.,{\bf 111}, 877946.\\



3418 



3419 
\noindent



3420 
{\bf [4]} Bluestein,\,H.\,B., 1993:



3421 
SynopticDynamic Meteorology in Midlatitudes. Volume II: Observations and Theory of Weather Systems, Oxford University Press, {594 pp}.\\



3422 



3423 
\noindent



3424 
{\bf [5]} Charney, J.\,B., 1993:



3425 
On the scale of atmospheric motions. Geofysisk.\,Publ.,{\bf 17}, No.\,2.\\



3426 



3427 
\noindent



3428 
{\bf [6]} Davis, C.\,A., 1992:



3429 
Piecewise potential vorticity inversion. J.\,Atmos.\,Sci., {\bf 49}, 13971411.\\



3430 



3431 
\noindent



3432 
{\bf [7]} Davis, C.\,A., and K.\,A.\,Emanuel, 1992:



3433 
Potential vorticity diagnostic of cyclogenesis. Mon.\,Wea.\,Rev., {\bf 119}, 19291953.\\



3434 



3435 
\noindent



3436 
{\bf [8]} Ertel, H., 1942:



3437 
Ein neuer hydrodynamischer Wirbelsatz. Meteor.\,Z., {\bf 59}, 277281.\\



3438 



3439 
\noindent



3440 
{\bf [9]} Fehlmann, R., 1997:



3441 
Dynamics of seminal PV elements.



3442 
PhD thesis No. 12229, ETH Z\"urich, 143 pp.\\



3443 



3444 
\noindent



3445 
{\bf [10]} Fehlmann, R. and H.\,C.\,Davies, 1997:



3446 
Misforecasts of synoptic systems: Diagnosis via PV retrodiction. Mon.\,Wea.\,Rev., {\bf 125}, 22472264.\\



3447 



3448 



3449 
\noindent



3450 
{\bf [11]} Holton, J.\,R., 1992:



3451 
An Introduction to Dynamic Meteorology. Third Edition.



3452 
Academic Press, San Diego, California, 511pp.\\



3453 



3454 
\noindent



3455 
{\bf [12]} Kleinschmidt, E., 1950:



3456 
\"Uber Aufbau und Entstehung von Zyklonen. Teil 1. Meteor.\,Rundsch., {\bf 3}, 16.\\



3457 



3458 
\noindent



3459 
{\bf [13]} Martius, O., E.\,Zenklusen, C. Schwierz, and H.\,C.\,Davies 2006:



3460 
Episodes of alpine heavy precipitation with an overlying elongated stratospheric intrusion: a climatology. Intl.\,J.\,Climatology., {\bf Vol. 96, Issue 9}, 11491164.\\



3461 



3462 
\noindent



3463 
{\bf [14]} McIntyre, M., 1997:



3464 
GEFD (Geophysical and Environmental Fluid Dynamics) Summer School, DAMPTP University of Cambridge.\\



3465 



3466 
\noindent



3467 
{\bf [15]} Lynch, P., 2006:



3468 
The Emergence of Numerical Weather Prediction. Richardson's Dream.



3469 
Cambridge UP, Cambridge, UK. 279pp.\\



3470 



3471 
\noindent



3472 
{\bf [16]} Nebeker, F., 1995:



3473 
Calculating the Weather. Meteorology in the 20th Century..



3474 
Academic Press, San Diego, California, 255pp.\\



3475 



3476 
\noindent



3477 
{\bf [17]} Press, W.\,H., Flannery, B.\,P., Teukolsky, S.\,A., and Vetterling, W.\,T., 1992:



3478 
Numerical Recipes in FORTRAN: The Art of Scientific Computing.



3479 
Cambridge University Press, Cambridge, UK, 1447pp.\\



3480 



3481 



3482 
\noindent



3483 
{\bf [18]} Sprenger, M.~and H.~Wernli, 2003:



3484 
A northern hemispheric climatology of crosstropopause exchange for the ERA15



3485 
time period (19791993).



3486 
{\sl J.~Geophys.~Res.,} {\bf 108}(D12), 8521, doi:10.1029/2002JD002636.\\



3487 



3488 
\noindent



3489 
{\bf [19]} Stohl, A., P. Bonasoni, P. Cristofanelli, W. Collins, J. Feichter,



3490 
A. Frank, C. Forster, E. Gerasopoulos, H. G\"uggeler, P. James,



3491 
T. Kentarchos, H. KrompKolb, B. Kr\"uger, C. Land, J. Meloen,



3492 
A. Papayannis, A. Priller, P. Seibert, M. Sprenger, G. J. Roelofs,



3493 
H. E. Scheel, C. Schnabel, P. Siegmund, L. Tobler, T. Trickl,



3494 
H. Wernli, V. Wirth, P. Zanis, and C. Zerefos: 2003.



3495 
Stratospheretroposphere exchange: A review, and what we have learned from STACCATO.



3496 
J.\,Geophys.\,Res.,{\bf Vol. 108}, No. D12, 8516, doi:10.1029/2002JD002490.\\



3497 



3498 



3499 
\noindent



3500 
{\bf [20]} Wallace, J.\,M. and P.\,V.\,Hobbs, R., 1977:



3501 
Atmospheric Science. An Introductory Survey.



3502 
Academic Press, San Diego, California, 467pp.\\



3503 



3504 



3505 



3506 
\newpage



3507 
\section{Appendix}



3508 



3509 
\subsection{Scalar coordinate transformation}



3510 



3511 
Four real function are needed for the transformation between the geographical latitude/longitude grid



3512 
used by ECMWF and the quasicartesian grid used by the PV inversion. Physically the new quasicartesian



3513 
coordinate system is obtained by introducing a new latitude/longitude grid, but now the two coordinates



3514 
being defined relative to a rotated north pole. Therefore, the new longitude and latitude are referred to as rotated coordinates. The position (in geographical latitude and longitude) of the new north pole is given as parameters POLPHI and POLLAM. Then, the functions\\[0.3cm]



3515 
\begin{center}



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



3517 
\begin{verbatim}



3518 
LAM = REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)



3519 
PHI = REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)



3520 
\end{verbatim}



3521 
\end{minipage}



3522 
\end{center}



3523 



3524 
\vspace{0.3cm}



3525 
\noindent



3526 
relate the rotated latitude (PHIS) and longitude (LAMS) to corresponding latitude (PHI) and



3527 
longitude (LAM) in the geographical system. Correspondingly, the following two functions give



3528 
the reverse transformation, i.e. to every position in geographical latitude/longitude a rotated



3529 
pair of coordinates is attributed. Again, the pole position (in the geographical grid) must be



3530 
passed as parameters (POLPHI and POLLAT).\\[0.3cm]



3531 
\begin{center}



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



3533 
\begin{verbatim}



3534 
LAMS = REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)



3535 
PHIS = REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)



3536 
\end{verbatim}



3537 
\end{minipage}



3538 
\end{center}



3539 



3540 
\vspace{0.3cm}



3541 
\noindent



3542 
Here, the geographical coordinates (PHI, LAM) are transformed into rotated coordinates (PHIS, LAMS).



3543 
If the corresponding pair of coordinates are given, any scalar field can easily be transformed. The following



3544 
Fortran functions are taken from the source code from the numerical weather prediction model of the German/Swiss weather service.



3545 



3546 



3547 
\begin{small}



3548 
\begin{verbatim}



3549 
C 



3550 
REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)



3551 
C 



3552 
C



3553 
C**** LMSTOLM  FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER



3554 
C**** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)



3555 
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT



3556 
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)



3557 
C** AUFRUF : LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)



3558 
C** ENTRIES : KEINE



3559 
C** ZWECK : BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER



3560 
C** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)



3561 
C** IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT



3562 
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)



3563 
C** VERSIONS



3564 
C** DATUM : 03.05.90



3565 
C**



3566 
C** EXTERNALS: KEINE



3567 
C** EINGABE



3568 
C** PARAMETER: PHIS REAL GEOGR. BREITE DES PUNKTES IM ROT.SYS.



3569 
C** LAMS REAL GEOGR. LAENGE DES PUNKTES IM ROT.SYS.



3570 
C** POLPHI REAL WAHRE GEOGR. BREITE DES NORDPOLS



3571 
C** POLLAM REAL WAHRE GEOGR. LAENGE DES NORDPOLS



3572 
C** AUSGABE



3573 
C** PARAMETER: WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION



3574 
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)



3575 
C**



3576 
C** COMMON



3577 
C** BLOECKE : KEINE



3578 
C**



3579 
C** FEHLERBE



3580 
C** HANDLUNG : KEINE



3581 
C** VERFASSER: D.MAJEWSKI



3582 



3583 
REAL LAMS,PHIS,POLPHI,POLLAM



3584 



3585 
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /



3586 



3587 
ZSINPOL = SIN(ZPIR18*POLPHI)



3588 
ZCOSPOL = COS(ZPIR18*POLPHI)



3589 
ZLAMPOL = ZPIR18*POLLAM



3590 
ZPHIS = ZPIR18*PHIS



3591 
ZLAMS = LAMS



3592 
IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS  360.0



3593 
ZLAMS = ZPIR18*ZLAMS



3594 



3595 
ZARG1 = SIN(ZLAMPOL)*( ZSINPOL*COS(ZLAMS)*COS(ZPHIS) +



3596 
1 ZCOSPOL* SIN(ZPHIS)) 



3597 
2 COS(ZLAMPOL)* SIN(ZLAMS)*COS(ZPHIS)



3598 
ZARG2 = COS(ZLAMPOL)*( ZSINPOL*COS(ZLAMS)*COS(ZPHIS) +



3599 
1 ZCOSPOL* SIN(ZPHIS)) +



3600 
2 SIN(ZLAMPOL)* SIN(ZLAMS)*COS(ZPHIS)



3601 
IF (ABS(ZARG2).LT.1.E30) THEN



3602 
IF (ABS(ZARG1).LT.1.E30) THEN



3603 
LMSTOLM = 0.0



3604 
ELSEIF (ZARG1.GT.0.) THEN



3605 
LMSTOLAM = 90.0



3606 
ELSE



3607 
LMSTOLAM = 90.0



3608 
ENDIF



3609 
ELSE



3610 
LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)



3611 
ENDIF



3612 



3613 
RETURN



3614 
END



3615 



3616 
C 



3617 
REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)



3618 
C 



3619 
C



3620 
C**** PHSTOPH  FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER



3621 
C**** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM



3622 
C**** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT



3623 
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)



3624 
C** AUFRUF : PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)



3625 
C** ENTRIES : KEINE



3626 
C** ZWECK : BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER



3627 
C** EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM



3628 
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT



3629 
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)



3630 
C** VERSIONS



3631 
C** DATUM : 03.05.90



3632 
C**



3633 
C** EXTERNALS: KEINE



3634 
C** EINGABE



3635 
C** PARAMETER: PHIS REAL GEOGR. BREITE DES PUNKTES IM ROT.SYS.



3636 
C** LAMS REAL GEOGR. LAENGE DES PUNKTES IM ROT.SYS.



3637 
C** POLPHI REAL WAHRE GEOGR. BREITE DES NORDPOLS



3638 
C** POLLAM REAL WAHRE GEOGR. LAENGE DES NORDPOLS



3639 
C** AUSGABE



3640 
C** PARAMETER: WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION



3641 
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)



3642 
C**



3643 
C** COMMON



3644 
C** BLOECKE : KEINE



3645 
C**



3646 
C** FEHLERBE



3647 
C** HANDLUNG : KEINE



3648 
C** VERFASSER: D.MAJEWSKI



3649 



3650 
REAL LAMS,PHIS,POLPHI,POLLAM



3651 



3652 
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /



3653 



3654 
SINPOL = SIN(ZPIR18*POLPHI)



3655 
COSPOL = COS(ZPIR18*POLPHI)



3656 
ZPHIS = ZPIR18*PHIS



3657 
ZLAMS = LAMS



3658 
IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS  360.0



3659 
ZLAMS = ZPIR18*ZLAMS



3660 
ARG = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)



3661 



3662 
PHSTOPH = ZRPI18*ASIN(ARG)



3663 



3664 
RETURN



3665 
END



3666 



3667 
C 



3668 
REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)



3669 
C 



3670 
C



3671 
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%



3672 
C



3673 
C**** LMTOLMS  FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM



3674 
C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)



3675 
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT



3676 
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)



3677 
C** AUFRUF : LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)



3678 
C** ENTRIES : KEINE



3679 
C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF



3680 
C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM



3681 
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT



3682 
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)



3683 
C** VERSIONS



3684 
C** DATUM : 03.05.90



3685 
C**



3686 
C** EXTERNALS: KEINE



3687 
C** EINGABE



3688 
C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM



3689 
C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM



3690 
C** POLPHI REAL GEOGR.BREITE DES NPOLS DES ROT. SYSTEMS



3691 
C** POLLAM REAL GEOGR.LAENGE DES NPOLS DES ROT. SYSTEMS



3692 
C** AUSGABE



3693 
C** PARAMETER: WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION



3694 
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)



3695 
C**



3696 
C** COMMON



3697 
C** BLOECKE : KEINE



3698 
C**



3699 
C** FEHLERBE



3700 
C** HANDLUNG : KEINE



3701 
C** VERFASSER: G. DE MORSIER



3702 



3703 
REAL LAM,PHI,POLPHI,POLLAM



3704 



3705 
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /



3706 



3707 
ZSINPOL = SIN(ZPIR18*POLPHI)



3708 
ZCOSPOL = COS(ZPIR18*POLPHI)



3709 
ZLAMPOL = ZPIR18*POLLAM



3710 
ZPHI = ZPIR18*PHI



3711 
ZLAM = LAM



3712 
IF(ZLAM.GT.180.0) ZLAM = ZLAM  360.0



3713 
ZLAM = ZPIR18*ZLAM



3714 



3715 
ZARG1 =  SIN(ZLAMZLAMPOL)*COS(ZPHI)



3716 
ZARG2 =  ZSINPOL*COS(ZPHI)*COS(ZLAMZLAMPOL)+ZCOSPOL*SIN(ZPHI)



3717 
IF (ABS(ZARG2).LT.1.E30) THEN



3718 
IF (ABS(ZARG1).LT.1.E30) THEN



3719 
LMTOLMS = 0.0



3720 
ELSEIF (ZARG1.GT.0.) THEN



3721 
LMTOLMS = 90.0



3722 
ELSE



3723 
LMTOLMS = 90.0



3724 
ENDIF



3725 
ELSE



3726 
LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)



3727 
ENDIF



3728 



3729 
RETURN



3730 
END



3731 



3732 
C 



3733 
REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)



3734 
C 



3735 
C



3736 
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%



3737 
C



3738 
C**** PHTOPHS  FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI



3739 
C**** AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)



3740 
C**** IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT



3741 
C**** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)



3742 
C** AUFRUF : PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)



3743 
C** ENTRIES : KEINE



3744 
C** ZWECK : UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF



3745 
C** EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM



3746 
C** ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT



3747 
C** DIE WAHREN KOORDINATEN (POLPHI, POLLAM)



3748 
C** VERSIONS



3749 
C** DATUM : 03.05.90



3750 
C**



3751 
C** EXTERNALS: KEINE



3752 
C** EINGABE



3753 
C** PARAMETER: PHI REAL BREITE DES PUNKTES IM GEOGR. SYSTEM



3754 
C** LAM REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM



3755 
C** POLPHI REAL GEOGR.BREITE DES NPOLS DES ROT. SYSTEMS



3756 
C** POLLAM REAL GEOGR.LAENGE DES NPOLS DES ROT. SYSTEMS



3757 
C** AUSGABE



3758 
C** PARAMETER: ROTIERTE BREITE PHIS ALS WERT DER FUNKTION



3759 
C** ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)



3760 
C**



3761 
C** COMMON



3762 
C** BLOECKE : KEINE



3763 
C**



3764 
C** FEHLERBE



3765 
C** HANDLUNG : KEINE



3766 
C** VERFASSER: G. DE MORSIER



3767 



3768 
REAL LAM,PHI,POLPHI,POLLAM



3769 



3770 
DATA ZRPI18 , ZPIR18 / 57.2957795 , 0.0174532925 /



3771 



3772 
ZSINPOL = SIN(ZPIR18*POLPHI)



3773 
ZCOSPOL = COS(ZPIR18*POLPHI)



3774 
ZLAMPOL = ZPIR18*POLLAM



3775 
ZPHI = ZPIR18*PHI



3776 
ZLAM = LAM



3777 
IF(ZLAM.GT.180.0) ZLAM = ZLAM  360.0



3778 
ZLAM = ZPIR18*ZLAM



3779 
ZARG = ZCOSPOL*COS(ZPHI)*COS(ZLAMZLAMPOL) + ZSINPOL*SIN(ZPHI)



3780 



3781 
PHTOPHS = ZRPI18*ASIN(ZARG)



3782 



3783 
RETURN



3784 
END



3785 
\end{verbatim}



3786 
\end{small}



3787 



3788 
\newpage



3789 
\subsection{Controlling Linux Shell Script}



3790 



3791 
In this appendix the controlling Linux Shell script of the PV inversion is reproduced. This script



3792 
offers a userfriendly interface to the inversion. It is essentially split into five



3793 
different section: (1) installation and parameter settings, (2) preparatory steps, (3)



3794 
iterative solution of the PV inversion problem, (4) postprocessing, and (5) diagnostic



3795 
tools.\\



3796 



3797 
\begin{small}



3798 
\begin{verbatim}



3799 
#!/bin/csh



3800 



3801 
# Master Linux script for PV inversion



3802 
# Michael Sprenger / Winter 2006,2007



3803 



3804 
# 



3805 
# Set some variables and paths



3806 
# 



3807 



3808 
# Handling of input parameters



3809 
set step="help"



3810 
if ( $#argv == 1 ) then



3811 
set step=$1



3812 
endif



3813 
if ( $#argv == 2 ) then



3814 
set step=$1



3815 
set file1=$2



3816 
endif



3817 
if ( $#argv == 3 ) then



3818 
set step=$1



3819 
set file1=$2



3820 
set file2=$3



3821 
endif



3822 



3823 
# Set base directory for programmes



3824 
set bdir=/home/sprenger/PV_Inversion_Tool/real/



3825 



3826 
# Set name of the parameter file, the coastline file and the sample files



3827 
set parafile="${PWD}/inversion.param"



3828 
set coastfile="${bdir}/prep/coastline.dat"



3829 
set sampledir="/net/rossby/lhome/sprenger/PV_Inversion_Tool/real/inp/"



3830 



3831 
# Extract parameters from parameter file



3832 
set progex=${bdir}/inversion.perl



3833 
set date=`${progex} ${parafile} date  awk '{ print $2}'`



3834 
set nofiter=`${progex} ${parafile} n_of_iteration  awk '{ print $2}'`



3835 
set save=`${progex} ${parafile} save_iteration  awk '{ print $2}'`



3836 
set idir=`${progex} ${parafile} inp_dir  awk '{ print $2}'`



3837 
set rdir=`${progex} ${parafile} run_dir  awk '{ print $2}'`



3838 
set odir=`${progex} ${parafile} out_dir  awk '{ print $2}'`



3839 



3840 
# 



3841 
# Installation, help and sample



3842 
# 



3843 



3844 
# Create the needed directories



3845 
if ( ${step} == "inst" ) then



3846 
if ( ! d ${idir} ) mkdir ${idir}



3847 
if ( ! d ${rdir} ) mkdir ${rdir}



3848 
if ( ! d ${odir} ) mkdir ${odir}



3849 
endif



3850 



3851 
# Create the needed directories



3852 
if ( ${step} == "help" ) then



3853 
echo



3854 
echo "Installation"



3855 
echo " inst: Creates the input, run and output directory"



3856 
echo



3857 
echo "Sample case study"



3858 
echo " sample: Copy all files for a sample case study"



3859 
echo



3860 
echo "Preparing input files [prep]"



3861 
echo " prep0: Calculate S file with PV and TH"



3862 
echo " prep1: Interpolate onto height levels"



3863 
echo " prep2: Rotate into local cartesian coordinate system"



3864 
echo " prep3: Add TH,PV,NSQ and RHO to the data file"



3865 
echo " prep4: Define modified and anomaly PV field and boundary values"



3866 
echo " prep5: Reduce the domain size and split the input files"



3867 
echo " prep6: Add the reference profile"



3868 
echo " prep7: Add coastlines to REF file"



3869 
echo " prep8: Move the files to the run directory"



3870 
echo



3871 
echo "Perform the PV inversion [pvin]"



3872 
echo " pvin1: Add NSQ, TH, RHO, and PV"



3873 
echo " pvin2: Change Ertel's PV anomaly into one of quasigeostrophic PV"



3874 
echo " pvin3: Inversion of quasigeostrophic PV anomaly"



3875 
echo " pvin4: Subtract anomaly from MOD file"



3876 
echo " pvin5: Keep iterative steps if save flag is set"



3877 
echo



3878 
echo "Postprocessing [post]"



3879 
echo " post1: Copy needed files from input and run directory"



3880 
echo " post2: Rotate from quasicartesian coordinate frame to lat/lon system"



3881 
echo " post3: Bring modified fields back to P file"



3882 
echo " post4: Calculate S file with PV and TH"



3883 
echo " post5: Make clean"



3884 
echo



3885 
echo "Diagnostic Tools"



3886 
echo " diag1: Check the consistency of the boundary conditions"



3887 
echo " diag2: Calculate the quasigeostrophic PV"



3888 
echo " diag3: Get the difference between two files"'



3889 
echo



3890 
endif



3891 



3892 



3893 
# Copy sample files (if specified)



3894 
if ( ${step} == "sample" ) then



3895 
\cp ${sampledir}/P20060116_18 ${idir}



3896 
\cp ${sampledir}/ml_cst ${idir}



3897 
\cp ${sampledir}/Z20060116_18 ${idir}



3898 
\cp ${sampledir}/pl_cst ${idir}



3899 
endif



3900 



3901 
# 



3902 
# Preparatory steps



3903 
# 



3904 



3905 
# Change to data directory



3906 
cd ${idir}



3907 



3908 
# Step 0: Calculate S file with PV and TH (not in prep mode)



3909 
if ( ${step} == "prep0" ) then



3910 
\rm f S${date}



3911 
p2s P${date} TH PV



3912 
endif



3913 



3914 
# Step 1: Interpolate onto height levels



3915 
if ( ${step} == "prep1"  ${step} == "prep" ) then



3916 
echo "P${date}" >! fort.10



3917 
echo "Z${date}" >> fort.10



3918 
echo "H${date}" >> fort.10



3919 
${bdir}/inversion.perl ${parafile} p2z >> fort.10



3920 
echo "U U P${date}" >> fort.10



3921 
echo "V V P${date}" >> fort.10



3922 
echo "T T P${date}" >> fort.10



3923 
echo "Q Q P${date}" >> fort.10



3924 
${bdir}/prep/p2z



3925 
\mv H${date}_cst zl_cst



3926 
changecst H${date} zl_cst



3927 
#\rm f fort.10



3928 
endif



3929 



3930 
# Step 2: Rotate into local cartesian coordinate system



3931 
if ( ${step} == "prep2"  ${step} == "prep" ) then



3932 
echo "H${date}" >! fort.10



3933 
echo "R${date}" >> fort.10



3934 
${bdir}/inversion.perl ${parafile} rotate_grid >> fort.10



3935 
echo "5" >> fort.10



3936 
echo "ORO" >> fort.10



3937 
echo "U.V" >> fort.10



3938 
echo "P" >> fort.10



3939 
echo "T" >> fort.10



3940 
echo "Q" >> fort.10



3941 
${bdir}/prep/rotate_grid



3942 
\mv R${date}_cst ro_cst



3943 
changecst R${date} ro_cst



3944 
\rm f fort.10



3945 
endif



3946 



3947 
# Step 3: Add TH,PV,NSQ and RHO to the data file



3948 
if ( ${step} == "prep3"  ${step} == "prep" ) then



3949 
echo "TH" >! fort.10



3950 
echo "R${date}" >> fort.10



3951 
echo "R${date}" >> fort.10



3952 
echo "5 " >> fort.10



3953 
${bdir}/prep/z2s



3954 
echo "PV" >! fort.10



3955 
echo "R${date}" >> fort.10



3956 
echo "R${date}" >> fort.10



3957 
echo "5 " >> fort.10



3958 
${bdir}/prep/z2s



3959 
echo "NSQ" >! fort.10



3960 
echo "R${date}" >> fort.10



3961 
echo "R${date}" >> fort.10



3962 
echo "5 " >> fort.10



3963 
${bdir}/prep/z2s



3964 
echo "RHO" >! fort.10



3965 
echo "R${date}" >> fort.10



3966 
echo "R${date}" >> fort.10



3967 
echo "5 " >> fort.10



3968 
${bdir}/prep/z2s



3969 
\rm f fort.10



3970 
endif



3971 



3972 
# Step 4: Set the modified PV field and boundary values



3973 
if ( ${step} == "prep4"  ${step} == "prep" ) then



3974 
echo "R${date}" >! fort.10



3975 
${bdir}/inversion.perl ${parafile} def_anomaly >> fort.10



3976 
${bdir}/prep/def_anomaly



3977 
\rm f fort.10



3978 
endif



3979 



3980 
# Step 5: Reduce the domain size and split the input files



3981 
if ( ${step} == "prep5"  ${step} == "prep" ) then



3982 
echo "PV PV R${date} ORG_${date}" >! fort.10



3983 
echo "U U R${date} ORG_${date}" >> fort.10



3984 
echo "V V R${date} ORG_${date}" >> fort.10



3985 
echo "TH TH R${date} ORG_${date}" >> fort.10



3986 
echo "Q Q R${date} ORG_${date}" >> fort.10



3987 
echo "P P R${date} ORG_${date}" >> fort.10



3988 
echo "T T R${date} ORG_${date}" >> fort.10



3989 
echo "PV_FILT PV_AIM R${date} MOD_${date}" >> fort.10



3990 
echo "U U R${date} MOD_${date}" >> fort.10



3991 
echo "V V R${date} MOD_${date}" >> fort.10



3992 
echo "Q Q R${date} MOD_${date}" >> fort.10



3993 
echo "P P R${date} MOD_${date}" >> fort.10



3994 
echo "T T R${date} MOD_${date}" >> fort.10



3995 
echo "NSQ NSQ R${date} MOD_${date}" >> fort.10



3996 
echo "RHO RHO R${date} MOD_${date}" >> fort.10



3997 
echo "TH TH R${date} MOD_${date}" >> fort.10



3998 
echo "PV_ANOM PV R${date} ANO_${date}" >> fort.10



3999 
echo "TH_ANOM TH R${date} ANO_${date}" >> fort.10



4000 
echo "UU_ANOM U R${date} ANO_${date}" >> fort.10



4001 
echo "VV_ANOM V R${date} ANO_${date}" >> fort.10



4002 
echo "ORO ORO R${date} REF_${date}" >> fort.10



4003 
echo "X X R${date} REF_${date}" >> fort.10



4004 
echo "Y Y R${date} REF_${date}" >> fort.10



4005 
echo "LAT LAT R${date} REF_${date}" >> fort.10



4006 
echo "LON LON R${date} REF_${date}" >> fort.10



4007 
echo "CORIOL CORIOL R${date} REF_${date}" >> fort.10



4008 
${bdir}/prep/cutnetcdf



4009 
\rm f fort.10



4010 
endif



4011 



4012 
# Step 6: Add the reference profile



4013 
if ( ${step} == "prep6"  ${step} == "prep" ) then



4014 
\rm f fort.10



4015 
echo "MOD_${date}" >! fort.10



4016 
echo "REF_${date}" >> fort.10



4017 
${bdir}/prep/ref_profile



4018 
\rm f fort.10



4019 
endif



4020 



4021 
# Step 7: Add coastlines to REF file



4022 
if ( ${step} == "prep7"  ${step} == "prep" ) then



4023 
\rm f fort.10



4024 
echo \"REF_${date}\" >! fort.10



4025 
echo \"${coastfile}\" >> fort.10



4026 
${bdir}/inversion.perl ${parafile} coastline >> fort.10



4027 
${bdir}/prep/coastline



4028 
\rm f fort.10



4029 
endif



4030 



4031 
# Step 8: Move the files to the run directory



4032 
if ( ${step} == "prep8"  ${step} == "prep" ) then



4033 
\mv MOD_${date} MOD_${date}_cst ${rdir}



4034 
\mv ORG_${date} ORG_${date}_cst ${rdir}



4035 
\mv ANO_${date} ANO_${date}_cst ${rdir}



4036 
\mv REF_${date} REF_${date}_cst ${rdir}



4037 
\rm f fort.10



4038 
endif



4039 



4040 
# 



4041 
# Inversion



4042 
# 



4043 



4044 
# Change to data directory



4045 
cd ${rdir}



4046 



4047 
# Start loop



4048 
set count=0



4049 
loop:



4050 



4051 
# Step 1: Add NSQ, TH, RHO, and PV to MOD file, take grid from REF file



4052 
if ( ${step} == "pvin1"  ${step} == "pvin" ) then



4053 
echo "NSQ" >! fort.10



4054 
echo "MOD_${date}" >> fort.10



4055 
echo "REF_${date}" >> fort.10



4056 
echo "5 " >> fort.10



4057 
${bdir}/pvin/z2s



4058 
echo "RHO" >! fort.10



4059 
echo "MOD_${date}" >> fort.10



4060 
echo "REF_${date}" >> fort.10



4061 
echo "5 " >> fort.10



4062 
${bdir}/pvin/z2s



4063 
echo "TH" >! fort.10



4064 
echo "MOD_${date}" >> fort.10



4065 
echo "REF_${date}" >> fort.10



4066 
echo "5 " >> fort.10



4067 
${bdir}/pvin/z2s



4068 
echo "PV" >! fort.10



4069 
echo "MOD_${date}" >> fort.10



4070 
echo "REF_${date}" >> fort.10



4071 
echo "5 " >> fort.10



4072 
${bdir}/pvin/z2s



4073 
\rm f fort.10



4074 
endif



4075 



4076 
# Step 2: Change Ertel's PV anomaly into an anomaly of quasigeostrophic PV



4077 
if ( ${step} == "pvin2"  ${step} == "pvin" ) then



4078 
echo "MOD_${date}" >! fort.10



4079 
echo "REF_${date}" >> fort.10



4080 
echo "ANO_${date}" >> fort.10



4081 
${bdir}/pvin/pv_to_qgpv



4082 
\rm f fort.10



4083 
endif



4084 



4085 
# Step 3: Inversion of quasigeostrophic PV anomaly with Neumann boundary



4086 
if ( ${step} == "pvin3"  ${step} == "pvin" ) then



4087 
echo "ANO_${date}" >! fort.10



4088 
echo "REF_${date}" >> fort.10



4089 
${bdir}/pvin/inv_cart



4090 
\rm f fort.10



4091 
endif



4092 



4093 
# Step 4: Prepare the output of the inversion for next iteration step



4094 
if ( ${step} == "pvin4"  ${step} == "pvin" ) then



4095 
\rm f fort.10



4096 
echo "MOD_${date}" >! fort.10



4097 
echo "ANO_${date}" >> fort.10



4098 
${bdir}/inversion.perl ${parafile} prep_iteration >> fort.10



4099 
${bdir}/pvin/prep_iteration



4100 
\rm f fort.10



4101 
endif



4102 



4103 
# Step 5: Keep iterative steps if save flag is set



4104 
if ( ${step} == "pvin5"  ${step} == "pvin" ) then



4105 
if ( "${save}" == "yes" ) then



4106 
set pre=''



4107 
if ( ${count} < 10 ) then



4108 
set pre='0'



4109 
endif



4110 
\cp MOD_${date} MOD_${date}_${pre}${count}



4111 
\cp ANO_${date} ANO_${date}_${pre}${count}



4112 
endif



4113 
endif



4114 



4115 
# End loop for iterations



4116 
if ( ${step} == "pvin" ) then



4117 
@ count = ${count} + 1



4118 
if ( ${count} < ${nofiter} ) goto loop



4119 
endif



4120 



4121 
# 



4122 
# Postprocessing



4123 
# 



4124 



4125 
# Change to data directory



4126 
cd ${odir}



4127 



4128 
# Step 1: Copy needed files from input and run directory



4129 
if ( ${step} == "post1"  ${step} == "post" ) then



4130 
ln sf ${rdir}/MOD_${date} ${rdir}/MOD_${date}_cst .



4131 
\cp ${idir}/P${date} ${idir}/ml_cst .



4132 
endif



4133 



4134 
# Step 2: Rotate from quasicartesian coordinate frame to lat/lon system



4135 
if ( ${step} == "post2"  ${step} == "post" ) then



4136 
echo "MOD_${date}" >! fort.10



4137 
echo "GEO_${date}" >> fort.10



4138 
${bdir}/inversion.perl ${parafile} rotate_lalo >> fort.10



4139 
echo "3" >> fort.10



4140 
echo "T" >> fort.10



4141 
echo "U.V" >> fort.10



4142 
echo "P" >> fort.10



4143 
${bdir}/post/rotate_lalo



4144 
\rm f fort.10



4145 
endif



4146 



4147 
# Step 3: Bring modified fields back to P file



4148 
if ( ${step} == "post3"  ${step} == "post" ) then



4149 
\rm f fort.10



4150 
echo "P${date}" >! fort.10



4151 
echo "GEO_${date}" >> fort.10



4152 
${bdir}/inversion.perl ${parafile} add2p >> fort.10



4153 
${bdir}/post/add2p



4154 
\rm f fort.10



4155 
endif



4156 



4157 
# Step 4: Calculate S file with PV and TH



4158 
if ( ${step} == "post4"  ${step} == "post" ) then



4159 
\rm f S${date}



4160 
p2s P${date} TH PV



4161 
endif



4162 



4163 
# Step 5: Make clean



4164 
if ( ${step} == "post5"  ${step} == "post" ) then



4165 
\rm f MOD_${date} MOD_${date}_cst



4166 
\rm f GEO_${date} GEO_${date}_cst



4167 
endif



4168 



4169 
# 



4170 
# Diagnostic Tools



4171 
# 



4172 



4173 
# Change to run directory



4174 
cd ${rdir}



4175 



4176 
# Step 1: Check the consistency of the boundary conditions (diag1)



4177 
if ( ${step} == "diag1" ) then



4178 
\rm f fort.10



4179 
echo "ANO_${date}" >! fort.10



4180 
echo "REF_${date}" >> fort.10



4181 
${bdir}/diag/check_boundcon



4182 
\rm f fort.10



4183 
endif



4184 



4185 
# Step 2: Calculate the quasigeostrophic PV (diag2 [ORGMODANO]



4186 
if ( ${step} == "diag2" ) then



4187 
\rm f fort.10



4188 
echo "${file1}_${date}" >! fort.10



4189 
echo "REF_${date}" >> fort.10



4190 
${bdir}/diag/calc_qgpv



4191 
\rm f fort.10



4192 
endif



4193 



4194 
# Step 3: Get difference between two files (diag3 [ORGMODANO]  [ORGMODANO])



4195 
if ( ${step} == "diag3" ) then



4196 
\rm f fort.10



4197 
echo "${file1}_${date}" >! fort.10



4198 
echo "${file2}_${date}" >> fort.10



4199 
echo "DIA_${date}" >> fort.10



4200 
${bdir}/diag/difference



4201 
\rm f fort.10



4202 
endif



4203 
\end{verbatim}



4204 
\end{small}



4205 



4206 
\newpage



4207 
\subsection{Controlling Parameter File}



4208 



4209 
The inversion problem is specified in a parameter file. This file is set as a parameter in the



4210 
controlling Linux Shell script (see Appendix~9.2) and is often the only file which has to be adapted.



4211 



4212 
\begin{small}



4213 
\begin{verbatim}



4214 
BEGIN DATA



4215 
DATE = 20060116_18; ! Date for case study



4216 
INP_DIR = /lhome/sprenger/PV_Inversion_Tool/real/inp; ! Input directory



4217 
RUN_DIR = /lhome/sprenger/PV_Inversion_Tool/real/run; ! Run directory



4218 
OUT_DIR = /lhome/sprenger/PV_Inversion_Tool/real/out; ! Output directory



4219 
END DATA



4220 



4221 
BEGIN GRID



4222 
GEO_XMIN = 180.; ! Geographical grid in zonal direction



4223 
GEO_NX = 361 ;



4224 
GEO_DX = 1.;



4225 
GEO_YMIN = 0.; ! Geographical grid in meridional direction



4226 
GEO_NY = 91 ;



4227 
GEO_DY = 1.;



4228 
GEO_ZMIN = 0.; ! Vertical levels ( ZMIN and DZ in [m] )



4229 
GEO_NZ = 125 ;



4230 
GEO_DZ = 200.;



4231 
ROT_NX = 250 ; ! Rotated grid ( DX and DY in [deg] )



4232 
ROT_NY = 250 ;



4233 
ROT_DX = 0.25;



4234 
ROT_DY = 0.25;



4235 
CLON = 65.; ! Longitude, latitude [deg] and angle of rotated grid



4236 
CLAT = 45.;



4237 
CROT = 0.;



4238 
END GRID



4239 



4240 
BEGIN ANOMALY



4241 
BOX_XMIN = 1200.; ! Box where to apply the digital filter



4242 
BOX_XMAX = 1200.;



4243 
BOX_YMIN = 1200.;



4244 
BOX_YMAX = 1200.;



4245 
BOX_
