diff --git a/Cat_UG_00/cat_ug.pdf b/Cat_UG_00/cat_ug.pdf
index e0b5abfa052065778bc25bae76ef8efa88bdc216..db104446bd4f258cd645d56e05395ec9d43e823c 100644
Binary files a/Cat_UG_00/cat_ug.pdf and b/Cat_UG_00/cat_ug.pdf differ
diff --git a/Cat_UG_00/chap_evoleq.tex b/Cat_UG_00/chap_evoleq.tex
index 12305becf9fc6d36d9413c2740efda781081a17d..567c07980ed8da9c118f4a3d5b96b56570cb2113 100644
--- a/Cat_UG_00/chap_evoleq.tex
+++ b/Cat_UG_00/chap_evoleq.tex
@@ -1,4 +1,4 @@
-%
+
 \chapter{Evolution equations of incompressible 2D fluids} 
 %
 %
@@ -20,7 +20,7 @@ From homogeneity and incompressibility of the fluid we get
   \rho \mathbf{\nabla \cdot u} = 0,
 \end{equation} 
 with $\rho$ the fluid density. Using this property we can introduce 
-a volume (mass) stream function measuring the volume (mass) flux  
+a volume (mass) streamfunction measuring the volume (mass) flux  
 across an arbitrary line from the point $(x_{0},y_{0})$ to a point $(x,y)$ 
 via the path integral
 \begin{equation} \label{eq_psizetadef}
@@ -45,10 +45,10 @@ via the path integral
 \end{equation}
 The minus sign in front of the integral just changes the direction
 of positive massflux across the line and is chosen to make 
-the classical stream function of 2D fluids compatible to 
-the stream function (see e.g.\ \cite{danilovandgurarie2000}) 
+the classical streamfunction of 2D fluids compatible to 
+the streamfunction (see e.g.\ \cite{danilovandgurarie2000}) 
 typically used for 2D rotating fluids in geosciences. 
-In terms of the stream function $\psi$ the integrated mass flux 
+In terms of the streamfunction $\psi$ the integrated mass flux 
 $\mathcal{M}$ across a line joining the points $(x_{0},y_{0})$ and 
 $(x,y)$ is given by
 \begin{equation} \label{eq_calM}
@@ -57,109 +57,235 @@ $(x,y)$ is given by
        \left[\psi(x,y) -  \psi(x_{0},y_{0}) \right],
 \end{equation}
 with $H_{0}$ the depth of the fluid. At the same time
-the stream function $\psi$ is connected to the velocity 
+the streamfunction $\psi$ is connected to the velocity 
 field $(u,v)$ and the (relative) 
 vorticity $\zeta = v_{x} - u_{y}$ via
 \begin{equation} \label{eq_psiuv}
   (u,v) = (- \psi_{y},\psi_{x}) \ \mbox{and} \ \ \zeta = \Delta \psi.
 \end{equation}
-Using the stream function $\psi$, the evolution equation (\ref{eq_2Dvortvel}
-can be also written in the form
+Using the streamfunction $\psi$, the evolution equation (\ref{eq_2Dvortvel})
+can be written in the form
 \begin{equation} \label{eq_2Dvortstream}
-  \zeta_{t} + J(\psi,\zeta) = F + D.
+  \zeta_{t} + J(\psi,\zeta) = F + D,
 \end{equation}
+with $J(\psi,\zeta) = \psi_{x} \zeta_{y} - \psi_{y} \zeta_{x}$ the Jacobian.
 Since the velocity field is divergence free (equation \ref{eq_conti})
-we can further derive the equivalent equation in flux form
+different representations of the Jacobian $J$ can be found. Using 
+appropriate superpositions of different representations one can derive 
+discrete analogs of the Jacobian $J$ which keep the symmetry 
+property $J(\psi,\zeta) = -J(\zeta,\psi)$ and conserve the integrated 
+vorticity $\zeta$, kinetic energy $\nabla \psi \cdot \nabla \psi $  and
+enstrophy $\zeta^{2}$. In \cite{arakawa1966} such a discrete analog is 
+derived for a finite difference numerical model. Keeping the notation of 
+\cite{arakawa1966} the Jacobian $J$ of equation (\ref{eq_2Dvortstream}) 
+is denoted by $J_{1}$. Starting from $J_{1}$ we can derive a second form 
+of the Jacobian 
+\begin{equation} \label{eq_Jacobi2} 
+  J_{2}(\psi,\zeta) 
+   = 
+  J_{1}(\psi,\zeta) + \psi \left( \zeta_{xy} - \zeta_{yx} \right)
+   = 
+  \left(\psi  \zeta_{y} \right)_{x}
+   - 
+  \left(\psi  \zeta_{x} \right)_{y},
+\end{equation}
+which gives the evolution equation
+\begin{equation} \label{eq_2DJacobi2} 
+  \zeta_{t} + \left(\psi  \zeta_{y} \right)_{x}  
+            - \left(\psi  \zeta_{x} \right)_{y} = F + D.
+\end{equation}
+Further we can derive a third form of Jacobian $J_{3}$ as follows 
+\begin{equation} \label{eq_Jacobi3}
+  J_{3}(\psi,\zeta) 
+   = 
+  J_{1}(\psi,\zeta) + \zeta \left( \psi_{xy} - \psi_{yx} \right)
+   = 
+  \left(\zeta  \psi_{x} \right)_{y}
+   - 
+  \left(\zeta  \psi_{y} \right)_{x},
+\end{equation}
+which inserted in the evolution equation yields the 2D fluid equations
+in flux form
 \begin{equation} \label{eq_2Dflux}
-  \zeta_{t} + \nabla \cdot \left(\mathbf{u} \zeta \right) = F + D.
+  \zeta_{t} + 
+  \left(\zeta  \psi_{x} \right)_{y}
+   - 
+  \left(\zeta  \psi_{y} \right)_{x},
+   = 
+  \zeta_{t} + \nabla \cdot \left(\mathbf{u} \zeta \right) 
+   = F + D.
 \end{equation}
-Starting from the flux form one can after repeated application of the
-continuity equation (\ref{eq_conti}) finally derive an equation where
-the Jacobian only depends on the two velocity fields $(u,v)$
-\begin{equation} \label{eq_2Duv}
-  \zeta_{t} 
-   + \partial_{x} \partial_{y} \left( v^{2} - u^{2} \right) 
-   + \left(\partial^{2}_{x} - \partial^{2}_{y} \right) uv
+It is possible to derive a fourth form of the Jacobian $J_{4}$
+which is a funciton of the field $v^{2} - u^{2})$ and 
+$uv$. We start again from the first Jacobian $J_{1}$
+\begin{equation} \label{eq_Jacobi4}
+  J_{4}(\psi,\zeta) 
+   = 
+  J_{1}(\psi,\zeta) + 
+   \left(2 \zeta - \psi_{x} \right) \ \left( \psi_{xy} - \psi_{yx} \right)
+   = 
+  \left(v^{2} - u^{2} \right)_{xy}  
+   + 
+  \left(uv \right)_{xx-yy},
+\end{equation}
+which leads to the evolution equation
+\begin{equation} \label{eq_2DJacobi4}
+  \zeta_{t}+ \left( v^{2} - u^{2} \right)_{xy} + \left( uv \right)_{xx - yy} 
    = 
   F + D.
 \end{equation}
-This form allows in the pseudo-spectral method (see \ref{sec_evolfourier}) 
-to reduce the number of Fourier transforms per time-step from 
-$3$, i.e.\ ($u,v,\zeta$) to $2$, i.e.\ ($u,v$).  
-
+All four forms given above are equivalent so that mathemtically one form is
+sufficient to catch all properties of the equations. Of course also other
+forms can be derived and it is also possible to superpose serveral forms.   
+The different forms get important if one tries to find discrete 
+representations of the Jacobian in the numerical scheme. Each form leads 
+to a discrete representation with different conservation and symmetry 
+properties. For more details see the chapters on the conservation 
+properties (\ref{sec_conprops}) and on the pseudo-spectral method 
+(\ref{sec_evolfourier}). 
 %
 \section{Quasi-two-dimensional rotating case} \label{sec_quasi2Dcase}
 %
 Starting from the shallow water equation on the $\beta$-plane one can
 derive (see e.g. \cite{danilovandgurarie2000}) an equation
 describing a rotating barotropic quasi-two-dimensional fluid which is a
-generalization of the 2D-equation (\ref{eq_2Dvortvel}). For small Rossby 
-numbers $\mathrm{Ro} = U/Lf$, with $U$ a typical horizontal fluid velocity
-at the length scale $L$ considered and $f$ the local coriolis paramter 
-we get the potential vorticity (PV) $q$ in quasi-geostrophic (QG) approximation 
+generalization of the 2D-equation (\ref{eq_2Dvortvel}). The shallow
+water potential vorticity is defined by 
+\begin{equation} \label{eq_qshallow}
+  q_{s} = \frac{\zeta + f}{H},
+\end{equation}
+where $H(x,y,t)$ is the fluid depth and $f$ the planetary vorticity. 
+Decomposing the fluid depth into a mean depth $H_{0}$, a constant 
+bottom topography B(x,y) and a time-dependent
+depth deviation $h(x,y,t)$ we can write $H(x,y,t) = H_{0} + h(x,y,t) - B(x,y)$. 
+Inserting the decomposition of the fluid depth into the equation 
+(\ref{eq_qshallow}) we can expand the potential vorticity
+\begin{equation} \label{eq_qshallowdecomp}
+  q_{s} 
+   =
+  \frac{\zeta + f}{H_{0} + h - B}
+   =
+  \frac{1}{H_{0}} \ \frac{\zeta + f}{1 + \Delta H / H_{0}}
+   =
+  \frac{1}{H_{0}} \ 
+  \left( \zeta + f \right) \left(1 - \frac{\Delta H}{H_{0}} + \dots \right),
+\end{equation}
+with $\Delta H = h - B$. Keeping only linear terms we finally get the barotopic
+vorticity
+\begin{equation} \label{eq_qbaro}
+  q = H_{0} q_{s} = \zeta - \frac{f}{H_{0}} \left(h - B\right) + f.
+\end{equation}
+Here we assume that the depth deviations $\Delta H = h - B$ are much smaller
+than the mean depth $H_{0}$. For small Rossby numbers $\mathrm{Ro} = U/Lf$, 
+with $U$ a typical horizontal velocity scale, $L$ a typical horizontal 
+length scale of the fluid motion and $f$ the local coriolis parameter 
+we can (see again e.g.\ \cite{danilovandgurarie2000}) introduce 
+the streamfunction 
+\begin{equation} \label{eq_psibaro}
+  \psi(x,y,t) = \frac{g}{f} \ h(x,y,t),    
+\end{equation}
+with $g$ the gravity. Using this streamfunction and expanding the 
+coriolis parameter linarly we can write the barotropic 
+quasi-geostrophic (QG) potential vorticity (PV) 
+on the $\beta$-plane (\ref{eq_qbaro}) as 
 \begin{equation} \label{eq_qdef}
-  q = \left( \nabla^{2}- \alpha^{2} \right) \psi + f,
+  q = \left( \nabla^{2}- \alpha^{2} \right) \psi 
+   + f_{0} + \beta y + \frac{f_{0}}{H_{0}} \ B.
+\end{equation}
+The linear expansion of the Coriolis parameter 
+$f(\varphi) = 2 \Omega \sin \varphi$ which is defined on a sphere 
+with radius $a$ and rotation rate $\Omega$ leads to
+\begin{equation} \label{eq_fbeta}
+  f(\varphi_{0} +  \Delta \varphi)
+   =
+  f(\varphi_{0}) + f^{\prime}(\varphi_{0}) \Delta \ \varphi   
+   = 
+  f_{0} + \beta  y,  
+\end{equation}
+with $\varphi_{0}$ the central latitude and the meridional 
+coordinate $y = a \Delta \varphi$ which is the linearization of 
+the projection $y = a \sin(\Delta \varphi)$. The $\beta$-parameter 
+is defined by $\beta = 2 \Omega \cos(\varphi_{0})/a$. The modification 
+parameter $\alpha = 1/L_{\mathrm{R}}$ is connected to the Rossby-Obukhov 
+radius of deformation $L_{\mathrm{R}} = \sqrt{g H_{0}}/f$. 
+As in the $2$-dimensional case the streamfunction $\psi$ 
+(remind the different definition) is related to the velocity field 
+$(u,v)$ again (see also equation \ref{eq_psiuv})  via
+\begin{equation} \label{eq_upsibaro}
+  u = -\psi_{y} \ \ \ \mbox{and} \ \ \ v = \psi_{x}.
 \end{equation}
-where we have the modification parameter $\alpha = 1/L^{2}_{\mathrm{R}}$ 
-with the Rossby-Obukhov radius of deformation 
-$L_{\mathrm{R}} = \sqrt{g H_{0}}/f$ and the stream function 
-$\psi = gh/f$. Here $g$ is the gravitational acceleration and $h(x,y)$ 
-the deviation of the mean fluid depth $H_{0}$ of the original 
-shallow water layer. As in the $2$-dimensional case the stream function $\psi$ 
-(remind the different definition) is related to the velocity fields $(u,v)$ 
-as defined in equation (\ref{eq_psiuv}). In an unforced non-dissipative fluid
-the QG PV is materially conserved
+In an unforced and non-dissipative fluid the QG PV is materially conserved
 \begin{equation} 
   q_{t} + J(\psi,q) = 0. 
 \end{equation}  
-Using the linear approximation of the coriolis parameter $f = f_{0} + \beta y$
+Using the linear approximation of the coriolis parameter (\ref{eq_fbeta})
 and introducing again forcing and dissipation we can write the evolution
-equation in the form
+equation for a barotropic fluid on the $\beta$-plane in the form
 \begin{equation} \label{eq_quasi2Dbaro}
-  q_{t} + J(\psi,q) + \beta \psi_{x} = F + D,
+  q_{t} + J(\psi,q + \frac{f_{0}}{H_{0}} \ B) + \beta \psi_{x} = F + D,
 \end{equation}
-with the vorticity $q$ given by
+with the vorticity $q$ again given by
 \begin{equation} \label{eq_vortquasi2Dbaro}
   q = \left(\nabla^2 -\alpha^{2} \right) \psi 
     = \zeta - \alpha^{2} \psi.
 \end{equation}
-Further one can introduce a variable mean fluid depth $H(x,y)$, 
-which in the simple case of a linear slope in $y$-direction leads 
-to a topographic $\beta$-effect (see e.g.\ \cite{vanheist1994}). 
+Here we used the property of the Jacobian $J(f,f) = 0$ for all fields
+$f(x,y)$ on the fluid domain.
 
+Idealizing the bottom topography to a linear slope in y-direction 
+$B(x,y) = B_{y} y$ the fluid experiences in addition to the ambient
+planetary $\beta$-effect a so called topographic $\beta$-effect 
+(see e.g.\ \cite{vanheist1994}) and the evolution equation simplifies to
+\begin{equation} \label{eq_qbarotopobeta} 
+  q_{t} + J(\psi,q) + \left(\beta + \frac{f_{0}}{H_{0}} B_{y} \right) \psi_{x} 
+   = F + D.
+\end{equation}
 In the form (\ref{eq_quasi2Dbaro}) and (\ref{eq_vortquasi2Dbaro}) 
-one can simulate incompressible 2D fluids and rotating quasi-2D fluids 
-with the same set of equations using different parameters. 
+one can simulate incompressible 2D fluids and rotating barotropic 
+quasi-2D fluids with the same set of equations using different parameters. 
 In this more general frame the simplest case of a non-rotating
 2D incompressible fluid is characterized by a vanishing 
 ambient vorticity gradient, i.e.\ $\beta = 0$, and the limit of
 an infinite Rossby radius $L_{\mathrm{R}} \longrightarrow \infty$
 or a vanishing modification parameter $\alpha \longrightarrow 0$.
+One has to keep in mind that the streamfunctions are different in 
+the two cases (see e.g.\ \cite{johnstonandliu2004}) and that there 
+are more subtle differences between 2D and QG Turbulence 
+(see e.g.\ \cite{tungandorlando2003}). 
+%
+\section{Multi-layer quasi-geostrophic case} \label{sec_multilayerqg}
+%
+%
+\section{The surface geostrophic case} \label{sec_sqg}
+%
+%
+\section{Conservation properties} \label{sec_conprops}
+%
+The properties of the Jacobian and the conditions at the fluid boundaries
+are at the base of the conservation properties of the fluid motions.
+Given two functions $g(x,y)$ and $h(x,y)$ defined on the fluid domains
+then the following integrals
+\begin{equation} \label{eq_intjacobian01}
+  \int_{x=0}^{X} \int_{x=0}^{Y} J(A,B) \ dxdy 
+  , \
+  \int_{x=0}^{X} \int_{x=0}^{Y} A J(A,B) \ dxdy 
+   \ \ \mbox{and} \ \    
+  \int_{x=0}^{X} \int_{x=0}^{Y} B J(A,B) \ dxdy 
+\end{equation}
+can be transformed through partial integration.
 
-One has to keep in mind that the stream functions are different in 
-the two cases. In the non-rotating case $\psi$ is defined by 
-equation (\ref{eq_calM}). In the rotating case we get 
-\begin{equation} \label{eq_hquasi2Dbaropsi}
-  h(x,y) = \frac{f}{g} \ \psi(x,y),    
-\end{equation} 
-so that $\psi$ is proportional to pressure deviations, which is not 
-the case in the non-rotating 2D case where the relation is more
-complex, see e.g.\ \cite{johnstonandliu2004}.
 
-Using the property of the Jacobian $J(f,f) = 0$ for all fields
-$f(x,y)$ on the fluid domain, equation (\ref{eq_quasi2Dbaro})
-is equivalent to
-\begin{equation} \label{eq_quasi2Dbarozeta}
-  q_{t} + J(\psi,\zeta) + \beta \psi_{x} = F + D,
-\end{equation}
-where the vorticity $q$ is still defined by equation 
-(\ref{eq_vortquasi2Dbaro}). From this from it directly follows that
-that the form of the 2D Jacobian in equations 
-(\ref{eq_2Dflux}) and (\ref{eq_2Duv}) can be also applied in the 
-quasi-two-dimensional rotating case.
 
-\section{Non-adiabatic terms}
+Above results can be generalized to more general fluid domains 
+(see e.g.\ \cite{salmonandtalley1989}).
+
+
+     
 
+
+%
+\section{Non-adiabatic terms}
+%
 \subsection{Laplacian based Viscosity and friction}
 
 Internal viscosity and external friction of the fluid are described by 
diff --git a/Cat_UG_00/chap_pseudospec.tex b/Cat_UG_00/chap_pseudospec.tex
index 330129dda1bf7a71169f35d37839aab187847f09..259dee39b709dc4c017d21b0fd71f616c99b1898 100644
--- a/Cat_UG_00/chap_pseudospec.tex
+++ b/Cat_UG_00/chap_pseudospec.tex
@@ -1079,11 +1079,11 @@ purely spectral and is called pseudo-spectral method, see e.g.\
 the higher wave numbers have to be filtered out. In CAT trunction is 
 used, see above. 
 
-We present three different forms of the Jacobian $J$ in physical space, 
-see equations (\ref{eq_2Dvortstream}), (\ref{eq_2Dflux}) 
-and (\ref{eq_2Duv}). Using the pseudo-spectral method for each
-form we get a different representation of the Jacobian $\hat{J}$ 
-in spectral space.   
+We present four different forms of the Jacobian $J$ in physical space, 
+see equations (\ref{eq_2Dvortstream}), (\ref{eq_2Dflux}), 
+(\ref{eq_2DJacobi2}) and (\ref{eq_2DJacobi4}). Using the pseudo-spectral 
+method for each form we get a different representation of the 
+Jacobian $\hat{J}$ in spectral space.   
 
 For the Jacobian (\ref{eq_2Dvortstream}) of the first form
 \begin{equation} \label{eq_jacobian01}
@@ -1091,18 +1091,18 @@ For the Jacobian (\ref{eq_2Dvortstream}) of the first form
    = 
   J(\psi,q) 
    = 
-  \partial_{x} \psi \ \partial_{y} q 
+  \psi_{x} \ q_{y} 
    -   
-  \partial_{x} q \ \partial_{y} \psi
+  q_{x} \  \psi_{y}
    =
-  \left[ v \ \partial_{y} q  + u \ \partial_{x} q \right],
+   v \ q_{y}  + u \  q_{x},
 \end{equation}
 one proceeds as follows. First the individual differential terms 
 are determined in spectral space using the fourier transform of 
 vorticity, and the spectral representation of the differential 
-operators. We get
+operators. The basic operators are given by 
 \begin{eqnarray} \label{eq_jacobian01_termsa}
-  \mathcal{F}(\partial_{x} \psi)
+  \mathcal{F}(\psi_{x})
   &=& 
   i  k_{x} \ \mathcal{F}(\psi)
    =
@@ -1111,11 +1111,11 @@ operators. We get
   - \ i \ \frac{k_{x}}{k^{2}_{x} + k^{2}_{y} + \alpha} \ 
   \mathcal{F}(q),
   \ \ 
-  \mathcal{F}(\partial_{y} q)
+  \mathcal{F}(q_{y})
    = 
   i k_{y} \ \mathcal{F}(q) \ \ \
   \\ \label{eq_jacobian01_termsb}
-  \mathcal{F}(\partial_{y} \psi)
+  \mathcal{F}(\psi_{y})
   &=&
   i  k_{y} \ \mathcal{F}(\psi)
    =
@@ -1124,7 +1124,7 @@ operators. We get
   - \ i \ \frac{k_{y}}{k^{2}_{x} + k^{2}_{y} + \alpha} \ 
   \mathcal{F}(q),
   \ \
-  \mathcal{F}(\partial_{x} q)
+  \mathcal{F}(q_{x})
    = 
   i k_{x} \ \mathcal{F}(q). \ \ \ 
 \end{eqnarray}
@@ -1137,43 +1137,41 @@ $\hat{J}$ is given by
 \begin{equation} \label{eq_jacobian01_J}
   \hat{J} 
    = 
-  \mathcal{F}(v \partial_{y} q)
-   -
-  \mathcal{F}(u \partial_{x} q).
+  \mathcal{F}(v q_{y})
+   +
+  \mathcal{F}(u q_{x})
+   = 
+  \mathcal{F}(v q_{y} + u q_{x}).
 \end{equation}
 Combining all necessary steps we can write
-\begin{eqnarray} \nonumber
+\begin{equation} \label{eq_jacobian01_Jall}
   \hat{J} 
-   &=&
+   =
   \mathcal{F}
-   \left(
+   \left[
     \mathcal{F}^{-1} 
      \left(
-      \ - \ i \frac{k_{x}}{k^{2}_{x}+k^{2}_{y}+\alpha} \ 
+       - i \frac{k_{x}}{k^{2}_{x}+k^{2}_{y}+\alpha} \ 
       \hat{q}
      \right) 
     \mathcal{F}^{-1} 
      \left(
       ik_{y} \ \hat{q}
      \right)
-   \right)
-    \\ \label{eq_jacobian01_Jall}
-   &-&
-  \mathcal{F}
-   \left(
+       +
     \mathcal{F}^{-1} 
      \left(
-     \ i \frac{k_{y}}{k^{2}_{x}+k^{2}_{y}+\alpha} \ 
+      i \frac{k_{y}}{k^{2}_{x}+k^{2}_{y}+\alpha} \ 
       \hat{q}
      \right) 
     \mathcal{F}^{-1}
      \left(
       ik_{x} \ \hat{q}
      \right)
-   \right),
-\end{eqnarray}
+   \right],
+\end{equation}
 where $\hat{q}$ is the vorticity in spectral space, the starting point
-for a new time step. As can be seen one needs $6$ $2$-FFT operations
+for a new time step. As can be seen one needs $5$ $2D$-FFT operations
 to determine the Jacobian. The components of the Jacobian 
 $\hat{J}_{k}$ are then used to determine the time 
 evolution of the different wave number components of the vorticity 
@@ -1182,11 +1180,7 @@ $\hat{q}_{\mathbf{k}}$, see equation (\ref{eq_evolqhat}).
 In the flux form of the evolution equation ({\ref{eq_2Dflux}})
 the second form of the Jacobian arises
 \begin{equation} \label{eq_jacobian02}
-  J(\psi,q)
-   = 
-  \partial_{x} (u \ q) 
-   +
-  \partial_{x} (v \ q).
+  J(\psi,q) = (u \ q)_{x} + (v \ q)_{y}.
 \end{equation}
 Following again equations (\ref{eq_jacobian01_termsa}) and 
 (\ref{eq_jacobian01_termsb}) we determine $\mathcal{F}(u)$ and
@@ -1235,14 +1229,47 @@ Combining again all necessary steps we can write
 As we can see the Jacobian in spectral space $\hat{J}$ can now be 
 determined by $5$ FFT operations.
 
-The third form of the Jacobian used in equation (\ref{eq_2Duv}) is
-given by
+The third form of the Jacobian in equation (\ref{eq_2DJacobi2})
+is given by
 \begin{equation} \label{eq_jacobian03}
+  J(\psi,q) = (\psi \ q_{y})_{x} - (\psi \ q_{x})_{y}.
+\end{equation}
+In spectral space we get 
+\begin{equation} \label{eq_jacobian03_J}
+  \hat{J} = -ik_{x} \ \mathcal{F}(\psi \ q_{y}) 
+            +ik_{y} \ \mathcal{F}(\psi \ q_{x}).
+\end{equation}
+Combining all necessary steps the Jacobian is determined by
+\begin{eqnarray} 
+  \hat{J} 
+   &=& 
+ -ik_{x} \ 
+  \mathcal{F}\left(
+     \mathcal{F}^{-1} 
+       \left(-\frac{1}{k_{x}^{2} + k_{y}^{2} + \alpha} \ \hat{q} \right) \
+     \mathcal{F}^{-1}
+       \left(
+         -ik_{y} \ \hat{q}
+       \right)
+             \right) 
+    \\ \label{eq_jacobian03_Jall}
+   &+&
+  ik_{y} \ 
+  \mathcal{F}\left(
+     \mathcal{F}^{-1} 
+       \left(-\frac{1}{k_{x}^{2} + k_{y}^{2} + \alpha} \ \hat{q} \right) \
+     \mathcal{F}^{-1}
+       \left(
+         -ik_{x} \ \hat{q}
+       \right)
+             \right).
+\end{eqnarray}
+The fourth form of the Jacobian used in equation (\ref{eq_2DJacobi4}) is \begin{equation} \label{eq_jacobian04}
   J(\psi,q)
    = 
   \partial_{x}\partial_{y} \left(v^{2} - u^{2} \right)
    +
-  \partial^{2}_{x} \partial^{2}_{y} uv.
+  \left( \partial^{2}_{x} - \partial^{2}_{y} \right)  uv.
 \end{equation}
 In this representation it is possible to reduce the number of FFT
 operations from $5$ to $4$. We first have to determine $\hat{u}$ 
@@ -1252,12 +1279,12 @@ space where the products $v^{2} - u^{2}$ and $uv$ are formed.
 Finally we have to transform them back to spectral space where they
 are differentiated. The Jacobian in spectral space $\hat{J}$ is 
 given by   
-\begin{equation} \label{eq_jacobian03_J}
+\begin{equation} \label{eq_jacobian04_J}
   \hat{J}
    = 
   -k_{x}k_{y} \ \mathcal{F}(v^{2} - u^{2})
    +
-   k^{2}_{x}k^{2}_{y} \mathcal{F}(uv)
+   \left( k^{2}_{x} - k^{2}_{y} \right) \mathcal{F}(uv)
 \end{equation} 
 or by combining all necessary steps
 \begin{eqnarray} \nonumber
@@ -1284,7 +1311,7 @@ or by combining all necessary steps
    \right)
     \\ \label{eq_jacobian02_Jall}
    &+&
-  k^{2}_{x} k^{2}_{y} \ 
+  \left( k^{2}_{x} - k^{2}_{y} \right) \ 
   \mathcal{F}
    \left(
     \mathcal{F}^{-1} 
diff --git a/Cat_UG_00/ref.bib b/Cat_UG_00/ref.bib
index 70185125b8b52045eff83b04b67552699d99f27d..76438b226f19a8e0c9d24c4bfa328debe57925b9 100644
--- a/Cat_UG_00/ref.bib
+++ b/Cat_UG_00/ref.bib
@@ -65,6 +65,14 @@
 @string{ pss =  {\it Planet.\ Space Sci.}}
 
 %==========================================================================
+@article{arakawa1966,
+   author    = {A. Arakawa},
+   title     = {{C}omputational {D}esign for {L}ong-{T}erm {N}umerical {I}ntegration of the {E}quations of {F}luid {M}otion: {T}wo-{D}imensional {I}ncompressible {F}low. {P}art {I}},
+   journal   = {J. Comput. Phys.},
+   pages     = {119--143 he},
+   year      = {1966},
+   volume    = {1}
+}
 
 @book{batchelor1967,
    author    = {G. K. Batchelor},
@@ -160,3 +168,20 @@
   pages =        {253--259}
 }
 
+@article{salmonandtalley1989,
+   author    = {R. Salmon and D. Talley},
+   title     = {Generalizations of Arakawa's Jacobian},
+   journal   = {\it Journal of Computational Physics},
+   year      = {1989},
+   volume    = {83}
+   pages     = {247--259}
+}
+
+@article{tungandorlando2003,
+   author    = {K. K. Tung and W. W. Orlando},
+   title     = {{O}n the Differences between 2D and QG Turbulence},
+   journal   = {\it Discrete and Continuous Dynamical Systems-Series B},
+   year      = {2003},
+   volume    = {3},
+   pages     = {145--162}
+}
diff --git a/cat/dat/GUI.cfg b/cat/dat/GUI.cfg
index 35d79c153b76bf83588635acfc5b5f31a8cdfbfb..91adfa92779aa18b2356f3b778ee44edc0757e21 100644
--- a/cat/dat/GUI.cfg
+++ b/cat/dat/GUI.cfg
@@ -9,21 +9,21 @@ Array:GQ
 Plot:ISOREC
 Palette:AUTO
 Title:Vorticity
-Geometry: 1000 1000    1    0
+Geometry: 1000 1000    1    1
 
 [Window 01]
-Array:GXM
-Plot:LINE
+Array:C4
+Plot:ISOAMP
 Palette:AUTO
-Title:GXM
-Geometry:  853  367  853   -1
+Title:Vorticity
+Geometry:  540  1000  1050   1
 
 [Window 02]
 Array:QC
 Plot:ISOAMP
 Palette:AMPLI
 Title:C4
-Geometry:  853  367 1706   -1
+Inactive:  853  367 1706   -1
 
 [Window 03]
 Array:GP
diff --git a/cat/dat/GUI_Amp.cfg b/cat/dat/GUI_Amp.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..91adfa92779aa18b2356f3b778ee44edc0757e21
--- /dev/null
+++ b/cat/dat/GUI_Amp.cfg
@@ -0,0 +1,139 @@
+Hamburg GUI Config File Version 16
+Screen: 2560x1418
+
+WinRows = 3
+WinCols = 3
+
+[Window 00]
+Array:GQ
+Plot:ISOREC
+Palette:AUTO
+Title:Vorticity
+Geometry: 1000 1000    1    1
+
+[Window 01]
+Array:C4
+Plot:ISOAMP
+Palette:AUTO
+Title:Vorticity
+Geometry:  540  1000  1050   1
+
+[Window 02]
+Array:QC
+Plot:ISOAMP
+Palette:AMPLI
+Title:C4
+Inactive:  853  367 1706   -1
+
+[Window 03]
+Array:GP
+Plot:MAPHOR
+Projection:POLAR
+Rotation factor:0
+Palette:P
+Title:Surface Pressure [hPa]
+Inactive:  853  367    0  389
+
+[Window 04]
+Array:GU
+Plot:MAPTRA
+Projection:AZIMUTHAL
+Rotation factor:0
+Palette:U
+Title:Tracer Level 3
+Inactive:  853  367  853  389
+
+[Window 05]
+Array:GV
+Plot:ISOLON
+Palette:U
+Title:Meridional Wind [m/s] Latitude 42N
+Inactive:  853  367 1706  389
+
+[Window 06]
+Array:GP
+Plot:ISOHOV
+Palette:P
+Title:Ps Hovmoeller Latitude 42N
+Inactive:  853  367    0  779
+
+[Window 07]
+Array:SCALAR
+Plot:ISOTS
+Palette:AUTO
+Title:Timeseries
+Inactive:  853  367  853  779
+
+[Window 08]
+Array:SCALAR
+Plot:ISOTAB
+Palette:AUTO
+Title:Tables
+Inactive:  853  367 1706  779
+
+[Control Window]
+Geometry:  920  122  820 1171
+
+# Scalar attributes for timeseries and table window
+
+[Scalar 00]
+Name:T
+Sub:
+Unit:C
+Scale:
+
+[Scalar 01]
+Name:RMS Div
+Sub:
+Unit:1/s
+Scale:-6
+
+[Scalar 02]
+Name:PE+KE
+Sub:
+Unit:m2/s2
+Scale:
+
+[Scalar 03]
+Name:RMS Ps
+Sub:
+Unit:hPa
+Scale:
+
+[Scalar 04]
+Name:Vort
+Sub:
+Unit:1/s
+Scale:
+
+# Parameter attributes for change menu
+
+[Parameter 00]
+ParName:NGUI
+ParInc:    1.0000
+ParMin:    1.0000
+ParMax: 1000.0000
+
+[Parameter 01]
+ParName:QMAX
+ParInc:    1.0000
+ParMin:    0.0000
+ParMax: 1000.0000
+
+[Parameter 02]
+ParName:DTNS
+ParInc:    5.0000
+ParMin:  -95.0000
+ParMax:  100.0000
+
+[Parameter 03]
+ParName:NSYNC
+ParInc:    1.0000
+ParMin:    0.0000
+ParMax:    1.0000
+
+[Parameter 04]
+ParName:EPSYNC
+ParInc:    1.0000
+ParMin:    0.0000
+ParMax:  100.0000
diff --git a/cat/src/cat.f90 b/cat/src/cat.f90
index 2a262e08a74a3b7dd5d642f67fa727f1312698a8..a429db51b775cf006f7dfc3723be81ef6b3f0be6 100644
--- a/cat/src/cat.f90
+++ b/cat/src/cat.f90
@@ -5,14 +5,14 @@ module catmod
 
 ! *********************************************
 ! *                                           *
-! *                 CAT                       * 
+! *                 CAT                       *
 ! *                                           *
 ! *        Computer Aided Turbulence          *
 ! *                                           *
 ! *        Version  1.1    July 2016          *
 ! *                                           *
 ! *********************************************
-! *  Theoretical Meteorology - KlimaCampus    * 
+! *  Theoretical Meteorology - KlimaCampus    *
 ! *         University of Hamburg             *
 ! *********************************************
 ! *     Hartmut Borth - Edilbert Kirk         *
@@ -42,13 +42,13 @@ module catmod
 character (256) :: catversion = "July 2016, Version 0.1"
 
 integer :: nshutdown = 0       ! Flag to stop program
-logical :: lrsttest  = .false. ! flag set to .true. for restart test, 
-                               ! parameter nsteps is not written 
+logical :: lrsttest  = .false. ! flag set to .true. for restart test,
+                               ! parameter nsteps is not written
                                ! to restart file. lrsttest is set to
                                ! true if a file RSTTEST is found in the
                                ! run directory. Call <make rsttest> to
                                ! do the restart test. To clean the run after
-                               ! a restart test call <make cleanresttest> 
+                               ! a restart test call <make cleanresttest>
 
 
 ! ***************************
@@ -109,20 +109,20 @@ integer, parameter :: ningp  = 5
 integer, parameter :: ninsp  = 5
 integer :: ingp(ningp) = &
    (/ qcde       , &
-      qfrccde    , &     
-       0         , &     
-       0         , &     
-       0           &    
-   /) 
+      qfrccde    , &
+       0         , &
+       0         , &
+       0           &
+   /)
 integer :: insp(ninsp) = &
    (/  qcde      , &
        qfrccde   , &
-       0         , &     
-       0         , &     
+       0         , &
+       0         , &
        0           &
-   /) 
+   /)
 
-!--- codes of variables to be written to cat_gp or cat_sp 
+!--- codes of variables to be written to cat_gp or cat_sp
 integer, parameter :: noutgp = 5
 integer, parameter :: noutsp = 5
 integer :: outgp(noutgp) = &
@@ -142,7 +142,7 @@ integer :: outsp(noutsp) = &
 
 
 ! *******************************
-! * Basic diagnostic parameters * 
+! * Basic diagnostic parameters *
 ! *******************************
 
 integer :: tsps          ! time steps per second (cpu-time)
@@ -164,7 +164,7 @@ real(8) :: ly = twopi   ! in y-direction (scale L_x = X/2pi)
 real(8) :: alpha = 0.0  ! alpha = 1/R_hat**2 [1/m**2]
                         ! with the non-dimensional Rossby
                         ! radius Ro_hat = Ro/L_x
-                            
+
 real(8) :: beta  = 0.0  ! ambient vorticity gradient [1/m*s]
 
 !---------------------!
@@ -176,20 +176,20 @@ integer :: diss_mthd  = 1  ! dissipation method
                            !    sig and lam and the powers psig and
                            !    plam
                            ! 2: Laplacian viscosity and friction
-                           !    characterized by the reciprocal 
+                           !    characterized by the reciprocal
                            !    damping time-scales rtsig and rtlam,
                            !    the cut-off wave-numbers ksig and klam
-                           !    and the powers psig and plam 
- 
+                           !    and the powers psig and plam
+
 
 !----------------------------------!
 ! Laplacian viscosity and friction !
 !----------------------------------!
 
-integer, parameter :: nsig = 2 ! maximum number of terms of Laplacian 
+integer, parameter :: nsig = 2 ! maximum number of terms of Laplacian
                                    ! dissipation on small scales
- 
-integer, parameter :: nlam = 2 ! maximum number of terms of Laplacian 
+
+integer, parameter :: nlam = 2 ! maximum number of terms of Laplacian
                                    ! dissipation on large scales
 
 !--- definition of dissipation on small scales
@@ -199,11 +199,11 @@ real(8) :: psig (nsig) =      & ! powers of Laplacian parameterizing
        /)
 real(8) :: sig(nsig)   =      & ! coefficients of different powers
        (/ 9.765625d-05,       & ! of the Laplacian [m^2*psig/s]
-          9.094947d-11        & 
+          9.094947d-11        &
        /)
 real(8) :: rtsig(nsig) =      & ! 1/time-scale of different powers
        (/ 1.d-1,              & ! of the Laplacian [1/s]
-          1.d+2               & 
+          1.d+2               &
        /)
 integer :: ksig (nsig) =      & ! lower "cut-off" wave number of different
        (/ 320,                & ! powers of Laplacian
@@ -214,14 +214,14 @@ integer :: ksig (nsig) =      & ! lower "cut-off" wave number of different
 real(8) :: plam (nlam) =      & ! powers of Laplacian modelling viscosity
        (/ -1.d0,              & ! on large scales plam <= 0
           -4.d0               &
-       /)      
+       /)
 real(8) :: lam(nlam)   =      & ! coefficients of different powers
        (/ 0.d0,               & ! of the Laplacian [m^(-2*psig)/s]
-          0.d0                & 
+          0.d0                &
        /)
 real(8) :: rtlam(nlam) =      & ! 1/time-scale of different powers
        (/ 0.d0,               & ! of the Laplacian [1/s]
-          0.d0                & 
+          0.d0                &
        /)
 integer :: klam (nlam) =      & ! lower "cut-off" wave number of different
        (/ 1,                  & ! powers of Laplacian
@@ -243,16 +243,16 @@ integer :: npert = 0 ! initial perturbation switch
                      ! 2 = white noise symmetric perturbation about the axis
                      !     y = pi (zero mean vorticity in the upper and
                      !     lower part of the fluid domain separately)
-                     ! 3 = white noise anti-symmetricy perturbation about 
-                     !     the axis y = pi (zero mean vorticity in the upper 
+                     ! 3 = white noise anti-symmetricy perturbation about
+                     !     the axis y = pi (zero mean vorticity in the upper
                      !     and lower power separately)
                      ! 4 = white noise symmetric perturbation about the axis
                      !     x = pi (zero mean vorticity in the upper and
                      !     lower part of the fluid domain separately)
-                     ! 5 = white noise anti-symmetricy perturbation about 
-                     !     the axis x = pi (zero mean vorticity in the upper 
+                     ! 5 = white noise anti-symmetricy perturbation about
+                     !     the axis x = pi (zero mean vorticity in the upper
                      !     and lower power separately)
- 
+
 real(8) :: apert = 1.0d-5 ! amplitude of perturbation
 
 
@@ -311,11 +311,11 @@ integer :: jacmthd  = 1       ! approximation method of Jacobian
                               ! 0 : no Jacobian
                               ! 1 : divergence form (uq)_x + (vq)_y
                               ! 2 : advection form   uq_x  + vq_y
-                              ! 3 : invariant form   
+                              ! 3 : invariant form
 
 !--- time integration
 integer :: nsteps    = 10000  ! number of time steps to be integrated
-integer :: tstep     = 0      ! current time step (since start of runs) 
+integer :: tstep     = 0      ! current time step (since start of runs)
 integer :: tstop     = 0      ! last time step of run
 
 integer :: tstp_mthd = 1      ! time stepping method
@@ -326,18 +326,18 @@ integer :: nstdout   = 2500   ! time steps between messages to
 
 
 integer :: ndiag     = 500    ! time steps between test outputs into
-                              ! out_diag (-1 no test outputs) 
+                              ! out_diag (-1 no test outputs)
 
 integer :: ngp       = 100     ! time steps between output of grid point
                                ! fields (-1 = no output)
 integer :: nsp       = 100     ! time steps between output of spectral
                                ! fields (-1 = no output)
-integer :: ncfl      = 100     ! time steps between cfl-check 
+integer :: ncfl      = 100     ! time steps between cfl-check
                                ! (-1 no cfl check)
 
 integer :: ntseri    = 100     ! time steps between time-series output
                                ! (-1 no time-series)
- 
+
 real(8) :: dt        = 1.d-3   ! length of time step [s]
 
 !--------------------------------------------!
@@ -358,7 +358,7 @@ real(8), allocatable :: guq(:,:)    ! u*q  [m/s^2]
 real(8), allocatable :: gvq(:,:)    ! v*q  [m/s^2]
 
 !--- case 2 (!! variables contain different fields depending on context)
-real(8), allocatable :: gqx (:,:) ! dq/dx [1/ms] 
+real(8), allocatable :: gqx (:,:) ! dq/dx [1/ms]
 real(8), allocatable :: gqy (:,:) ! dq/dy [1/ms]
 real(8), allocatable :: gjac(:,:) ! jacobian [s^-2]
 
@@ -372,7 +372,7 @@ real(4), allocatable :: ggui(:,:) ! single precision transfer
 
 !----------------------------------------------------------!
 ! variables in spectral space, representation as imaginary !
-! and real part (prefix f) used for fourier transformation ! 
+! and real part (prefix f) used for fourier transformation !
 ! and complex representation (prefix c) used for time      !
 ! evolution                                                !
 !----------------------------------------------------------!
@@ -424,10 +424,10 @@ complex(8), allocatable :: cli(:,:)    ! linear time propagation
 ! *******************
 
 !--- gui communication (guimod)
-integer :: ngui    = 1   ! global switch and parameter 
-                         ! ngui > 0 GUI = on  
-                         ! ngui specifies moreover the number of 
-                         ! time-steps between two calls of gui_transfer 
+integer :: ngui    = 1   ! global switch and parameter
+                         ! ngui > 0 GUI = on
+                         ! ngui specifies moreover the number of
+                         ! time-steps between two calls of gui_transfer
                          ! which exchanges data between GUI and CAT
 
 integer :: nguidbg = 0   ! GUI debug mode
@@ -437,9 +437,9 @@ real(4) :: parc(5) = 0.0 ! timeseries display
 
 
 !--- code of predefined simulations (simmod)
-integer :: nsim  = 0 ! 0 no predefined simulation is specified 
+integer :: nsim  = 0 ! 0 no predefined simulation is specified
                      ! nsim > 0 predefined simulation is run
-                     ! 
+                     !
                      ! available predefined simulations (experiments):
                      ! -----------------------------------------------
                      !    code       name
@@ -447,19 +447,19 @@ integer :: nsim  = 0 ! 0 no predefined simulation is specified
                      !     1    decaying_jet01
                      !    51    top-hat wind-forcing
                      !    52    top-hat wind-forcing
-                     ! ----------------------------------------------- 
-                     ! 
+                     ! -----------------------------------------------
+                     !
                      ! Predefined simulations are given in the dat
-                     ! directory as sim_XXXX.nl with XXXX the code  
+                     ! directory as sim_XXXX.nl with XXXX the code
                      ! To activate a given simulation one has to copy
-                     ! a given sim_XXXX.nl to sim_namelist in the run 
+                     ! a given sim_XXXX.nl to sim_namelist in the run
                      ! directory.
-                     ! ----------------------------------------------- 
+                     ! -----------------------------------------------
 
 !--- usermod (usermod)
 integer :: nuser = 0 ! 1/0 user mode is switched on/off.
-                     ! In user mode it is possible to introduce 
-                     ! new code into CAT (see chapter <Modifying 
+                     ! In user mode it is possible to introduce
+                     ! new code into CAT (see chapter <Modifying
                      ! CAT> in the User's Guide).
 
 
@@ -480,10 +480,10 @@ real(8), allocatable :: guym(:)   ! mean x-velocity allong y-dir [m/s]
 real(8), allocatable :: gvxm(:)   ! mean y-velocity allong x-dir [m/s]
 real(8), allocatable :: gvym(:)   ! mean y-velocity allong y-dir [m/s]
 
-real(4), allocatable :: gguixm(:) ! single precision mean in x-dir for gui 
-real(4), allocatable :: gguiym(:) ! single precision mean in y-dir for gui 
+real(4), allocatable :: gguixm(:) ! single precision mean in x-dir for gui
+real(4), allocatable :: gguiym(:) ! single precision mean in y-dir for gui
 
-real(4), allocatable :: ggxm(:,:) ! single precision mean in y-dir for gui 
+real(4), allocatable :: ggxm(:,:) ! single precision mean in y-dir for gui
 
 end module catmod
 
@@ -684,7 +684,7 @@ if (lcatnl) then
   read(nucatnl,cat_nl)
 endif
 
-!--- check consistency of namelist parameters   
+!--- check consistency of namelist parameters
 if (jacmthd .lt. 0 .or. jacmthd .gt. 3) then
    write(nudiag, &
    '(/," ************************************************")')
@@ -790,21 +790,21 @@ allocate(gpsi(1:ngx,1:ngy)) ; gpsi(:,:) = 0.0  ! stream function
 
 allocate(ggui(1:ngx,1:ngy)) ; ggui(:,:)  = 0.0  ! GUI transfer
 
-!--- spetral space 
-allocate(fu(0:nfx,0:nfy))   ; fu(:,:)   = 0.0 ! u 
-allocate(fv(0:nfx,0:nfy))   ; fv(:,:)   = 0.0 ! v 
+!--- spetral space
+allocate(fu(0:nfx,0:nfy))   ; fu(:,:)   = 0.0 ! u
+allocate(fv(0:nfx,0:nfy))   ; fv(:,:)   = 0.0 ! v
 
 allocate(k2(0:nkx,0:nfy))      ; k2(:,:)      = 0.0 ! Laplacian
 allocate(k2pa(0:nkx,0:nfy))    ; k2pa(:,:)    = 0.0 ! modified Laplacian
 allocate(rk2pa(0:nkx,0:nfy))   ; rk2pa(:,:)   = 0.0 ! modified Laplacian^-1
 allocate(kxrk2pa(0:nkx,0:nfy)) ; kxrk2pa(:,:) = 0.0 ! q --> v
 allocate(kyrk2pa(0:nkx,0:nfy)) ; kyrk2pa(:,:) = 0.0 ! q --> u
-allocate(kx2mky2(0:nkx,0:nfy)) ; kx2mky2(:,:)  = 0.0 ! utility for Jacobian 3 
-allocate(kxky  (0:nkx,0:nfy))  ; kxky  (:,:)  = 0.0 ! utility for Jacobian 3 
+allocate(kx2mky2(0:nkx,0:nfy)) ; kx2mky2(:,:)  = 0.0 ! utility for Jacobian 3
+allocate(kxky  (0:nkx,0:nfy))  ; kxky  (:,:)  = 0.0 ! utility for Jacobian 3
 
 
 
-allocate(cli(0:nkx,0:nfy))  ; cli(:,:)   = (0.0,0.0) ! linear time propagator 
+allocate(cli(0:nkx,0:nfy))  ; cli(:,:)   = (0.0,0.0) ! linear time propagator
 
 allocate(cq(0:nkx,0:nfy))   ; cq(:,:)    = (0.0,0.0) ! vorticity
 allocate(c4(0:nkx,0:nfy))   ; c4(:,:)    = (0.0,0.0) ! vorticity
@@ -838,7 +838,7 @@ case(1)
    allocate(guq(1:ngx,1:ngy)) ; guq(:,:) = 0.0  ! u*q
    allocate(gvq(1:ngx,1:ngy)) ; gvq(:,:) = 0.0  ! v*q
 
-   !--- spectral fields  
+   !--- spectral fields
    allocate(fuq(0:nfx,0:nfy)) ; fuq(:,:) = 0.0  ! u*q
    allocate(fvq(0:nfx,0:nfy)) ; fvq(:,:) = 0.0  ! v*q
 case(2)
@@ -849,15 +849,15 @@ case(2)
    allocate(gqy(1:ngx,1:ngy))  ; gqy(:,:)  = 0.0  ! qy
    allocate(gjac(1:ngx,1:ngy)) ; gjac(:,:) = 0.0  ! jacobian in grid-space
 
-   !--- spectral fields  
+   !--- spectral fields
    allocate(fqx(0:nfx,0:nfy)) ; fqx(:,:) = 0.0  ! u*qx or qx
    allocate(fqy(0:nfx,0:nfy)) ; fqy(:,:) = 0.0  ! v*qy or qy
 case(3)
    !--- allocate fields used by jacobian
-   !--- grid point fields 
-   allocate(gv2mu2(1:ngx,1:ngy)) ; gv2mu2(:,:) = 0.0  ! v^2-u^2 
-   allocate(guv   (1:ngx,1:ngy)) ; guv   (:,:) = 0.0  ! u*v 
-   !--- spetral fields  
+   !--- grid point fields
+   allocate(gv2mu2(1:ngx,1:ngy)) ; gv2mu2(:,:) = 0.0  ! v^2-u^2
+   allocate(guv   (1:ngx,1:ngy)) ; guv   (:,:) = 0.0  ! u*v
+   !--- spetral fields
    allocate(fv2mu2(0:nfx,0:nfy)) ; fv2mu2(:,:) = 0.0  ! v^2-u^2
    allocate(fuv   (0:nfx,0:nfy)) ; fuv   (:,:) = 0.0  ! u*v
    print *, "init_jacobian case 3"
@@ -879,34 +879,34 @@ select case (jacmthd)
 case (1)
    !--- deallocate fields used by jacobian
 
-   !--- grid point fields 
+   !--- grid point fields
    deallocate(guq)
    deallocate(gvq)
 
-   !--- spetral fields 
+   !--- spetral fields
    deallocate(fuq)
    deallocate(fvq)
 
 case (2)
    !--- deallocate fields used by jacobian
 
-   !--- grid point fields 
+   !--- grid point fields
    deallocate(gqx)
    deallocate(gqy)
    deallocate(gjac)
 
-   !--- spetral fields 
+   !--- spetral fields
    deallocate(fqx)
    deallocate(fqy)
 
 case (3)
    !--- deallocate fields used by jacobian
 
-   !--- grid point fields 
+   !--- grid point fields
    deallocate(gv2mu2)
    deallocate(guv)
 
-   !--- spetral fields 
+   !--- spetral fields
    deallocate(fv2mu2)
    deallocate(fuv)
 
@@ -1162,7 +1162,7 @@ integer    :: jx,jy,n
 !
 ! arg = -dt*[diss_sig + diss_lam + beta_term]
 !
-!  with 
+!  with
 !     diss_sig  = -sum_j=1,nsig sig(j)*k^(2*psig(j))
 !     diss_lam  = -sum_j=1,nlam lam(j)*k^(2*plam(j))
 !     b_term    = i*beta*kx/(k^2 + alpha)
@@ -1230,7 +1230,7 @@ integer :: jx,jy
 
 if (nforc .eq. 0) return
 
-!---     
+!---
 ent=0.d0
 do jy = 0, nfy
   do jx = 0, nkx
@@ -1239,10 +1239,10 @@ do jy = 0, nfy
     if (jy.ne.0.or.jx.ne.0) then
       if (sqrt(k2tmp).ge.kfmin.and.sqrt(k2tmp).le.kfmax) then
         nk = nk + 1
-        in(nk,1) = jx 
+        in(nk,1) = jx
         in(nk,2) = jy
 !       ent=ent+k2+alpha
-        ent=ent+k2pa(jy,jx)
+        ent=ent+k2pa(jx,jy)
       endif
     endif
   enddo
@@ -1281,7 +1281,7 @@ if (nforc .eq. 0) return
 
 !--- forcing is specified by input files
 
-!--- examples are produced in simmod.f90 
+!--- examples are produced in simmod.f90
 
 !--- produce masks
 
@@ -1328,7 +1328,7 @@ select case(npert)
       do jy = 1,ngy/2
          call random_number(xtmp)
          xtmp = 2.0*apert*(xtmp - sum(xtmp)*rnx)
-         gqpert(:,jy)       = xtmp(:)   
+         gqpert(:,jy)       = xtmp(:)
          gqpert(:,ngy+1-jy) = xtmp(:)
       enddo
       call grid_to_fourier(gqpert,cqpert,nfx,nfy,ngx,ngy)
@@ -1339,35 +1339,35 @@ select case(npert)
       do jy = 1,ngy/2
          call random_number(xtmp)
          xtmp = 2.0*apert*(xtmp - sum(xtmp)*rnx)
-         gqpert(:,jy)       =  xtmp(:)   
+         gqpert(:,jy)       =  xtmp(:)
          gqpert(:,ngy+1-jy) = -xtmp(:)
       enddo
       call grid_to_fourier(gqpert,cqpert,nfx,nfy,ngx,ngy)
       cqpert(0,0) = (0.0,0.0)
       cq = cq + cqpert
    case(4)
-      !--- symmetric about x = pi 
+      !--- symmetric about x = pi
       do jx = 1,ngx/2
          call random_number(ytmp)
          ytmp = 2.0*apert*(ytmp - sum(ytmp)*rny)
-         gqpert(jx,:)       = ytmp(:)   
+         gqpert(jx,:)       = ytmp(:)
          gqpert(ngx+1-jx,:) = ytmp(:)
       enddo
       call grid_to_fourier(gqpert,cqpert,nfx,nfy,ngx,ngy)
       cqpert(0,0) = (0.0,0.0)
       cq = cq + cqpert
    case(5)
-      !--- anti-symmetric about x = pi 
+      !--- anti-symmetric about x = pi
 !     do jx = 1,ngx/2
 !        call random_number(ytmp)
 !        ytmp = 2.0*apert*(ytmp - sum(ytmp)*rny)
-!        gqpert(jx,:)       =  ytmp(:)   
+!        gqpert(jx,:)       =  ytmp(:)
 !        gqpert(ngx+1-jx,:) = -ytmp(:)
 !     enddo
       do jy = 1,ngy/2
          call random_number(xtmp)
          xtmp = 2.0*apert*(xtmp - sum(xtmp)*rnx)
-         gqpert(:,jy)       =  xtmp(:)   
+         gqpert(:,jy)       =  xtmp(:)
          gqpert(:,ngy+1-jy) = -xtmp(:)
       enddo
       gqpert = transpose(gqpert)
@@ -1448,7 +1448,7 @@ if((ngp .gt. 0) .and. mod(tstep,ngp).eq.0)  then
       case (qcde)
          call cat_wrtgp(nugp,gq,qcde,0)
       end select
-   enddo 
+   enddo
 endif
 
 if((nsp .gt. 0) .and. mod(tstep,nsp).eq.0)  then
@@ -1480,13 +1480,18 @@ implicit none
 ggui(:,:) = gq(:,:) ! double precision -> single
 call guiput("GQ" // char(0), ggui, ngx, ngy, 1)
 
+! double -> single and center about ky = 0
+c4(:,0:nky-1)     = cq(:,nky+1:nfy)  
+c4(:,nky:nfy) = cq(:,0:nky)
+call guiput("C4" // char(0), c4, nkx+1, nfy+1, 2)   ! fc
+
 if (npost > 0) then
    !--- zonal means
    gguixm(:) = gqxm(:)   ! double -> single
-   call guiput("GQXM" // char(0), gguixm, ngy, 1, 1)   ! q 
+   call guiput("GQXM" // char(0), gguixm, ngy, 1, 1)   ! q
 
    gguixm(:) = gpsixm(:) ! double -> single
-   call guiput("GPSIXM" // char(0), gguixm, ngy, 1, 1) ! psi 
+   call guiput("GPSIXM" // char(0), gguixm, ngy, 1, 1) ! psi
 
 !  gguixm(:) = guxm(:) ! double -> single
 !  call guiput("GUXM" // char(0), gguixm, ngy, 1, 1)   ! u
@@ -1502,10 +1507,10 @@ if (npost > 0) then
 
    !--- meridional means
    gguiym(:) = gqym(:)   ! double -> single
-   call guiput("GQYM" // char(0), gguiym, ngx, 1, 1) ! q 
+   call guiput("GQYM" // char(0), gguiym, ngx, 1, 1) ! q
 
    gguiym(:) = gpsiym(:) ! double -> single
-   call guiput("GPSIYM" // char(0), gguiym, ngx, 1, 1) ! psi 
+   call guiput("GPSIYM" // char(0), gguiym, ngx, 1, 1) ! psi
 
    gguiym(:) = guym(:) ! double -> single
    call guiput("GUYM" // char(0), gguiym, ngx, 1, 1)   ! u
@@ -1513,9 +1518,6 @@ if (npost > 0) then
    gguiym(:) = gvym(:) ! double -> single
    call guiput("GVYM" // char(0), gguiym, 1, ngx, 1)   ! v
 
-   c4(:,:) = cq(:,:)
-   call guiput("C4" // char(0), c4, nkx+1, nfy+1, 2)   ! v
-
 endif
 
 return
@@ -1603,16 +1605,16 @@ subroutine init_tstep
 use catmod
 implicit none
 
-real(8) :: dt2 
+real(8) :: dt2
 
 
 dt2 = dt/2
 
 !--- determine tstop and return if tstep > 0
 if (tstep .eq. 0) then
-   tstop = tstep + nsteps      
+   tstop = tstep + nsteps
 else
-   tstop = tstep - 1 + nsteps   
+   tstop = tstep - 1 + nsteps
    return
 endif
 
@@ -1671,10 +1673,10 @@ allocate(guym(1:ngx))   ; guym(:)   = 0.0 ! u-mean in y-dir [m/s]
 allocate(gvxm(1:ngy))   ; gvxm(:)   = 0.0 ! v-mean in x-dir [m/s]
 allocate(gvym(1:ngx))   ; gvym(:)   = 0.0 ! v-mean in y-dir [m/s]
 
-allocate(gguixm(1:ngy)) ; gguixm(:) = 0.0 ! single precision mean in x-dir  
-allocate(gguiym(1:ngx)) ; gguiym(:) = 0.0 ! single precision mean in y-dir  
+allocate(gguixm(1:ngy)) ; gguixm(:) = 0.0 ! single precision mean in x-dir
+allocate(gguiym(1:ngx)) ; gguiym(:) = 0.0 ! single precision mean in y-dir
 
-allocate(ggxm(1:ngx,4)) ; ggxm(:,:) = 0.0 ! single precision mean in x-dir  
+allocate(ggxm(1:ngx,4)) ; ggxm(:,:) = 0.0 ! single precision mean in x-dir
 
 return
 end subroutine init_post
@@ -1688,11 +1690,11 @@ use catmod
 implicit none
 
 if (npost > 0) then
-   deallocate(gqxm) 
-   deallocate(gqym) 
+   deallocate(gqxm)
+   deallocate(gqym)
 
    deallocate(gpsixm)
-   deallocate(gpsiym) 
+   deallocate(gpsiym)
 
    deallocate(guxm)
    deallocate(guym)
@@ -1860,7 +1862,7 @@ call put_restart_iarray('npert',npert,1,1)
 call put_restart_array ('apert',apert,1,1)
 call put_restart_iarray('nforc',nforc,1,1)
 call put_restart_iarray('kfmin',kfmin,1,1)
-call put_restart_iarray('kfmax',kfmax,1,1) 
+call put_restart_iarray('kfmax',kfmax,1,1)
 call put_restart_array ('aforc',aforc,1,1)
 call put_restart_array ('tforc',tforc,1,1)
 call put_restart_iarray('myseed',myseed,mxseedlen,1)
@@ -1894,7 +1896,7 @@ subroutine close_files
 use catmod
 
 if (lrst)            close(nurstini)
-if (lcatnl)          close(nucatnl) 
+if (lcatnl)          close(nucatnl)
 close(nurstfin)
 close(nudiag)
 if (ntseri .gt. 0)   close(nutseri)
@@ -1930,7 +1932,7 @@ case (1)
 
    !--- (!!! Jacobian has different sign than in CAT-UG !!!)
    do jx = 0, nkx
-      cjac0(jx,:) = & 
+      cjac0(jx,:) = &
        cmplx( kx(jx)*fuq(jx+jx+1,:)+ky(:)*fvq(jx+jx+1,:), &
              -kx(jx)*fuq(jx+jx  ,:)-ky(:)*fvq(jx+jx  ,:))
    enddo
@@ -1943,18 +1945,18 @@ case (2)
 
 case (3)
    !--- Jacobian in invariant form (third representation in CAT-UG)
-   gv2mu2 = gv**2 - gu**2 
+   gv2mu2 = gv**2 - gu**2
    guv    = gu*gv
 
    call grid_to_fourier(gv2mu2,fv2mu2,nfx,nfy,ngx,ngy)
    call grid_to_fourier(guv,fuv,nfx,nfy,ngx,ngy)
 
    do jx = 0, nkx
-      cjac0(jx,:) = & 
+      cjac0(jx,:) = &
        cmplx(kxky(jx,:)*fv2mu2(2*jx  ,:) + kx2mky2(jx,:)*fuv(2*jx  ,:), &
              kxky(jx,:)*fv2mu2(2*jx+1,:) + kx2mky2(jx,:)*fuv(2*jx+1,:))
    enddo
-      
+
 end select
 
 if (jac_scale .ne. 1.0 .or. jac_scale .ne. 0.0)  &
@@ -2021,9 +2023,9 @@ eni=0.d0
 enf=0.d0
 
 select case (nforc)
-   case (1) 
+   case (1)
       cq = cq + cqfrc
- 
+
    case (2)
       do ifk = 1,nk
          i = in(ifk,1)
@@ -2126,7 +2128,7 @@ real(8) :: cfl
 real(8) :: maxu,maxv
 
 
-!--- check courant number 
+!--- check courant number
 maxu = maxval(abs(gu(:,:))) * dt * ngx / lx
 maxv = maxval(abs(gv(:,:))) * dt * ngy / ly
 
@@ -2145,10 +2147,10 @@ subroutine init_ops
 use catmod
 implicit none
 
-integer :: jx,jy    
- 
+integer :: jx,jy
+
 kx2    = kx*kx
-kx2(0) = 1.0    
+kx2(0) = 1.0
 
 ky2    = ky*ky
 
@@ -2157,7 +2159,7 @@ do jy = 0, nfy
       k2     (jx,jy) = (kx2(jx)+ky2(jy))        ! - Laplacian      psi --> -q
       k2pa   (jx,jy) = (kx2(jx)+ky2(jy))+alpha  ! - mod Laplacian  psi --> -q
       rk2pa  (jx,jy) = 1.0d0/(k2(jx,jy)+alpha)  ! 1/mod Laplacian    q --> -psi
-      kxrk2pa(jx,jy) = kx(jx)*rk2pa(jx,jy)      ! q ---> v/i  
+      kxrk2pa(jx,jy) = kx(jx)*rk2pa(jx,jy)      ! q ---> v/i
       kyrk2pa(jx,jy) = ky(jy)*rk2pa(jx,jy)      ! q ---> u/i
       kx2mky2(jx,jy) = kx2(jx)-ky2(jy)          ! utility Jacobian 3
       kxky   (jx,jy) = kx(jx)*ky(jy)            ! utility Jacobian 3
diff --git a/cat/src/guix11.c b/cat/src/guix11.c
index e2e33a5b5f60135126d17efad57c20302002032f..7f0de9d1497ed2c6bf7913f5220f34689eb4fd9b 100644
--- a/cat/src/guix11.c
+++ b/cat/src/guix11.c
@@ -1,5 +1,5 @@
 /*
-    guix11 - a GUI for PUMA, PLASIM & CAT 
+    guix11 - a GUI for PUMA, PLASIM & CAT
     2016 Edilbert Kirk & Hartmut Borth
 
     This program is free software; you can redistribute it and/or modify
@@ -925,7 +925,7 @@ void AzimuthalImage(struct MapImageStruct *s, struct MapImageStruct *d)
    int dy ;           // centre relative y position
    int p00;           // euqator
    int l00;           // reference longitude
-   
+
    unsigned int dih;  // destination image height
    unsigned int diw;  // destination image width
    unsigned int dpw;  // destination image padded width
@@ -1008,15 +1008,15 @@ void SwapIEEE16(char W[2])
 {
    char B;
 
-   B = W[0]; W[0] = W[1]; W[1] = B; 
+   B = W[0]; W[0] = W[1]; W[1] = B;
 }
 
 void SwapIEEE32(char W[4])
 {
    char B;
 
-   B = W[0]; W[0] = W[3]; W[3] = B; 
-   B = W[1]; W[1] = W[2]; W[2] = B; 
+   B = W[0]; W[0] = W[3]; W[3] = B;
+   B = W[1]; W[1] = W[2]; W[2] = B;
 }
 
 int ReadINT(FILE *fpi)
@@ -1102,8 +1102,8 @@ void ReadImage(struct MapImageStruct *ei, char *filename)
    BuffImageData = malloc(PicBytes);
    fseek(fp,OffBits,SEEK_SET);
    n = fread(BuffImageData,PicBytes,1,fp);
-   fclose(fp);  
-   
+   fclose(fp);
+
    ei->X = XCreateImage(display,CopyFromParent,ScreenD,ZPixmap,0,
                   ei->d,PadWidth,ImageBMI.Height,8,0);
 
@@ -1252,7 +1252,7 @@ int AllocateColorCells(struct ColorStrip cs[])
       XAllocNamedColor(display,colormap,cs[i].Name,&xcolor1,&xcolor2);
       cs[i].pixel = xcolor1.pixel;
       ++i;
-   }   
+   }
    return i;
 }
 
@@ -1373,13 +1373,13 @@ void TestWindow(int xx, int yy, int wi, int he, int *x, int *y,
    unsigned int DepthReturn  = 0;
    unsigned int BorderReturn = 0;
    Window TW,Parent,Child;
-   
+
    TW = XCreateWindow(display,RootWindow(display,screen_num),
         xx,yy,wi,he,                   // x,y,w,h
         0,CopyFromParent,InputOutput,  // border, depth, class
           CopyFromParent,0,NULL);      // visual, valuemask, attributes
    XSelectInput(display,TW,ExposureMask);
-   XMapWindow(display,TW);           
+   XMapWindow(display,TW);
    XMoveWindow(display,TW,xx,yy);
    XNextEvent(display,&CowEvent);            // Wait until WM mapped it
    XGetGeometry(display,TW,&Parent,x,y,W,H,&BorderReturn,&DepthReturn);
@@ -1398,7 +1398,7 @@ void CreateTestWindow(void)
 {
    int X,Y,xp,yp;
    unsigned int W,H;
-   
+
    // Open a test window with displacement x = 100 and y = 100
    // Test the returned coordinates for upper left corner
    // and compute left margin (border) and top margin (title bar)
@@ -1606,7 +1606,7 @@ void CreateControlWindow(void)
    strcpy(Title1,"Unknown model");
    if (Model == 0) // PUMA
    {
-      if (MRpid < 0) 
+      if (MRpid < 0)
          strcpy(Title1,"PUMA - KlimaCampus Hamburg");
       else
          sprintf(Title1,"Run %d: PUMA - KlimaCampus Hamburg",MRpid);
@@ -1836,7 +1836,7 @@ void ShowGridStatus(void)
    if (Grid) CheckMark(x,y,d);
 }
 
-      
+
 int RedrawControlWindow(void)
 {
    char Text[80];
@@ -1847,7 +1847,7 @@ int RedrawControlWindow(void)
    status = XGetWindowAttributes(display,Cow,&CurAtt);
    WinXSize = CurAtt.width;
    WinYSize = CurAtt.height;
-   
+
    XSetForeground(display,gc,BlackPix);
    XSetBackground(display,gc,WhitePix);
 
@@ -1865,7 +1865,7 @@ int RedrawControlWindow(void)
    for (i=0 ; i < PSDIM ; ++i)
    if (pixelstar[j][i] == '*')
       XDrawPoint(display,Cow,gc,i+x1,j+y1);
-  
+
    XSetForeground(display,gc,LightRed.pixel);
    x1 = 10; x2 = 40 ; y1 = OffsetY + 10 ; y2 = OffsetY + 40;
    XDrawLine(display,Cow,gc,x2  ,y1-1,x2  ,y2  );
@@ -1895,7 +1895,7 @@ int RedrawControlWindow(void)
    XSetForeground(display,gc,Blue.pixel);
    XFillPolygon(display,Cow,gc,PauseButton1,4,Convex,CoordModeOrigin);
    XFillPolygon(display,Cow,gc,PauseButton2,4,Convex,CoordModeOrigin);
-   
+
    x1 = 95; x2 = 102; y1 = OffsetY + 10; y2 = OffsetY + 40;
    XSetForeground(display,gc,DarkBlue.pixel);
    XDrawLine(display,Cow,gc,x1-1,y1-1,x1-1,y2  );
@@ -1959,7 +1959,7 @@ int CheckEndianess(void)
       int  i;
    } ec;
 
-   ec.i = 8; 
+   ec.i = 8;
    return (ec.b[0] == 0);
 }
 
@@ -1983,7 +1983,7 @@ void initgui_(int *model, int *debug, int *lats, int *mrpid, int *mrnum, char *p
 
    gettimeofday(&TimeVal,NULL);
    Seed = TimeVal.tv_sec;
-   
+
    Model = *model;
    Debug = *debug;
    Latitudes = *lats;
@@ -2000,7 +2000,7 @@ void initgui_(int *model, int *debug, int *lats, int *mrpid, int *mrnum, char *p
    if (Debug) printf("initgui(%d,%d)\n",Model,Debug);
    if ((display=XOpenDisplay(display_name)) == NULL)
    {
-      fprintf(stderr,"%s: cannot connect to X server %s\n", 
+      fprintf(stderr,"%s: cannot connect to X server %s\n",
               progname, XDisplayName(display_name));
       exit(1);
    }
@@ -2063,7 +2063,7 @@ void initgui_(int *model, int *debug, int *lats, int *mrpid, int *mrnum, char *p
    wm_hints.initial_state = NormalState;
    wm_hints.input = True;
    wm_hints.flags = StateHint | InputHint;
-   
+
    class_hints.res_name = progname;
    class_hints.res_class = "PUMA";
 
@@ -2115,7 +2115,7 @@ void initgui_(int *model, int *debug, int *lats, int *mrpid, int *mrnum, char *p
 
    for (i=0 ; i < NUMPAL ; ++i)
       LineCo[i] = AllocateColorCells(Pallet[i]);
-  
+
    /* Color cells for control window */
 
    XAllocNamedColor(display,colormap,"red"        ,&Red       ,&Dummy);
@@ -2646,7 +2646,7 @@ void LinePlot(int w)
    }
 
    // fill plot area with black
-   
+
    XSetBackground(display,gc,BlackPix);
    XSetForeground(display,gc,BlackPix);
    XFillRectangle(display,pix,gc,0,0,WinXSize,WinYSize);
@@ -2670,12 +2670,12 @@ void LinePlot(int w)
    for (j=0 ; j < DimY ; ++j)
    {
       // scale plot data
-   
+
       for (i=0 ; i < DimX ; ++i)
       {
          LIxp[w][i].y = InYSize - f * (Field[i+j*DimX] - zmin);
       }
-   
+
       XSetForeground(display,gc,Simplestrip[j].pixel);
       XDrawLines(display,pix,gc,LIxp[w],DimX,CoordModeOrigin);
    }
@@ -2952,7 +2952,7 @@ void TracerPlot(int w)
    {
       for (j=0 ; j < DimY; ++j)
       {
-         SpeedScale[j] = InXSize * DeltaTime * 60.0 / 40000000.0 
+         SpeedScale[j] = InXSize * DeltaTime * 60.0 / 40000000.0
                        / cos((j-0.5*(DimY-1)) * (M_PI/DimY));
       }
       SpeedScale[DimY-1] = SpeedScale[0] = SpeedScale[1];
@@ -3111,11 +3111,11 @@ void SH_Amplitudes(int w)
 }
 
 
-/* ========== */
-/* Amplitudes */
-/* ========== */
+/* ============= */
+/* FC_Amplitudes */
+/* ============= */
 
-void Amplitudes(int w)
+void FC_Amplitudes(int w)
 {
    int i,j,m,n;
    int x,y;
@@ -3155,14 +3155,17 @@ void Amplitudes(int w)
 
    dx = WinXSize / (DimX+1);
    dy = WinYSize / (DimY+1);
+   if (dx < dy) dy = dx;
+   else         dx = dy;
    r  = (dx+dx+dy+dy) / 6;
    XSetForeground(display,gc,BlackPix);
-   for (m=0,i=0 ; m < DimX ; ++m)
+   i = 0;
+   for (n=0 ; n < DimY ; ++n)
    {
-      y = FixFontHeight + 4 + m * dx;
-      for (n=0 ; n < DimY ; ++n)
+      y = FixFontHeight + 4 + n * dy;
+      for (m=0 ; m < DimX ; ++m)
       {
-         x = dx/2 + n * dx;
+         x = dx/2 + m * dx;
          XSetForeground(display,gc,Aco[w][i++].pixel);
          XFillArc(display,pix,gc,x,y,r,r,0,360*64);
       }
@@ -3201,7 +3204,7 @@ void ReopenWindow(int i)
    {
       if (WinAtt[i].x < ScreenOffset) WinAtt[i].x = ScreenOffset + OutXSize * (i % WinCols);
       if (WinAtt[i].y < ScreenOffset) WinAtt[i].y = ScreenOffset + OutYSize * (i / WinCols);
-      Win[i] = XCreateSimpleWindow(display,RootWindow(display,screen_num), 
+      Win[i] = XCreateSimpleWindow(display,RootWindow(display,screen_num),
                WinAtt[i].x,WinAtt[i].y,WinAtt[i].w,WinAtt[i].h,
                BorderWidth,BlackPix,WhitePix);
       XSetWMProtocols(display,Win[i],&Delwin,1);
@@ -3265,7 +3268,7 @@ void DisplayHelpWindow(void)
    h = BigFontHeight * (HelpLines + 1);
    x = (ScreenW - w) / 2;
    y = (ScreenH - h) / 2;
-   HelpWindow = XCreateSimpleWindow(display,RootWindow(display,screen_num), 
+   HelpWindow = XCreateSimpleWindow(display,RootWindow(display,screen_num),
                 x,y,w,h,BorderWidth,Yellow.pixel,DarkBlue.pixel);
    XStringListToTextProperty(&HelpTitle,1,&HelpName);
    XSetWMProtocols(display,HelpWindow,&Delwin,1);
@@ -3619,7 +3622,7 @@ void guiclose_(void)
    }
    while (!(XCheckTypedWindowEvent(display,Cow,ButtonPress,&CowEvent) &&
             HitStopButton(&CowEvent)));
-     
+
 // XCloseDisplay(display); // segmentation fault on sun compiler!
 }
 
@@ -3728,7 +3731,7 @@ void lp2ps(REAL *a, REAL *b)
    {
       dx = arx * (xnp - i);
       dy = ary * (ynp - j);
-      x = atan2(dx,dy) * lfa;      // angle 
+      x = atan2(dx,dy) * lfa;      // angle
       y = sqrt(dx * dx + dy * dy); // distance from NP
       if (x < 0.0) x += (DimX-1);
       k = x;                       // integer part
@@ -4028,23 +4031,23 @@ void ShowGridPolar(void)
       XSetForeground(display,gc,WhitePix);
       XSetBackground(display,gc,BlackPix);
       XSetLineAttributes(display,gc,1,LineOnOffDash,CapButt,JoinRound);
-   
+
       /* Northern Hemisphere */
-   
+
       XDrawArc(display,pix,gc,OffX+dx/6,OffY+dy/6,2*dx/3,2*dy/3,0,360*64);
       XDrawArc(display,pix,gc,OffX+dx/3,OffY+dy/3,dx/3,dy/3,0,360*64);
-   
+
       XDrawLine(display,pix,gc,OffX,yh,InXSize,yh);
       XDrawLine(display,pix,gc,OffX+xh,OffY,OffX+xh,InYSize);
-   
+
       /* Southern Hemisphere */
-   
+
       XDrawArc(display,pix,gc,OffX+7*dx/6,OffY+dy/6,2*dx/3,2*dy/3,0,360*64);
       XDrawArc(display,pix,gc,OffX+4*dx/3,OffY+dy/3,dx/3,dy/3,0,360*64);
-   
+
       x = OffX + dx + xh;
       XDrawLine(display,pix,gc,x,OffY,x,InYSize);
-   
+
       if (GridLabel)
       {
          x  = OffX + xh - FixFontWidth/2;
@@ -4206,7 +4209,7 @@ void iso(int w,int PicType,REAL *field,int dimx,int dimy,int dimz,int pal)
    char Text[128];
 
    int i,j,k,len,lens,xp,yp,status,x;
-   int y,dx,r,width,height;
+   int y,dx,dy,r,width,height;
    INTXU border,depth;
    REAL f,o,ra,rb;
    int CapLines;
@@ -4221,7 +4224,7 @@ void iso(int w,int PicType,REAL *field,int dimx,int dimy,int dimz,int pal)
    // At high frame rates (SkipFreq > 1) some types are redrawn
    // at "SkipFreq" intervals
 
-   if (SkipFreq > 1 && (nstep % SkipFreq) != 0 && 
+   if (SkipFreq > 1 && (nstep % SkipFreq) != 0 &&
       (PicType == ISOCS  || PicType == ISOHOR ||
        PicType == MAPHOR || PicType == ISOCOL ||
        PicType == ISOREC )) return;
@@ -4436,12 +4439,12 @@ void iso(int w,int PicType,REAL *field,int dimx,int dimy,int dimz,int pal)
    pix = WinPixMap[w].Pix; /* Set current pixmap */
 
    /* Draw colour bar */
-   
+
    if ((SizeChanged || Cstrip == Autostrip) &&
        (PicType == ISOCS  || PicType == ISOHOR ||
-       	PicType == ISOHOV || PicType == ISOLON ||
+        PicType == ISOHOV || PicType == ISOLON ||
         PicType == ISOCOL || PicType == MAPHOR ||
-        PicType == ISOREC ))
+        PicType == ISOREC || PicType == ISOAMP))
    {
       XSetForeground(display,gc,BlackPix);
       XFillRectangle(display,pix,gc,0,InYSize,WinXSize,WinYSize);
@@ -4489,7 +4492,7 @@ void iso(int w,int PicType,REAL *field,int dimx,int dimy,int dimz,int pal)
    }
 
    /* Draw mode legend */
-   
+
    if (SizeChanged && PicType == ISOSH)
    {
       dx = WinXSize / 23;
@@ -4552,6 +4555,73 @@ void iso(int w,int PicType,REAL *field,int dimx,int dimy,int dimz,int pal)
       }
    }
 
+   /* Double Fourier Coefficients */
+
+   if (SizeChanged && PicType == ISOAMP)
+   {
+      dx = WinXSize / (DimX+1);
+      dy = WinYSize / (DimY+1);
+      if (dx < dy) dy = dx;
+      else         dx = dy;
+      XSetForeground(display,gc,BlackPix);
+      XFillRectangle(display,pix,gc,0,0,WinXSize,WinYSize);
+
+      XSetForeground(display,gc,LightGreen.pixel);
+      XSetBackground(display,gc,BlackPix);
+      for (i=0 ; i < DimX ; i+=2)
+      {
+         sprintf(Text,"%d",i);
+         len    = strlen(Text);
+         width  = XTextWidth(FixFont,Text,len);
+         height = FixFont->ascent + FixFont->descent;
+         xp     = dx + i * dx - width/2 - FixFontWidth/2 + 1;
+         yp     = FixFontHeight;
+         XDrawImageString(display,pix,gc,xp,yp,Text,len);
+      }
+      XSetForeground(display,gc,LightBlue.pixel);
+      for (i=0 ; i < DimY ; i+=2)
+      {
+         sprintf(Text,"%d",i);
+         len    = strlen(Text);
+         width  = XTextWidth(FixFont,Text,len);
+         height = FixFont->ascent + FixFont->descent;
+         xp     = WinXSize - width - 2;
+         yp     = 2 * FixFontHeight + i * dy;
+         if (yp+FixFont->descent > InYSize) break;
+         XDrawImageString(display,pix,gc,xp,yp,Text,len);
+      }
+      strcpy(Text,"n/m");
+      len    = strlen(Text);
+      width  = XTextWidth(FixFont,Text,len);
+      xp     = WinXSize - width - 2;
+      yp     = FixFontHeight;
+      XSetForeground(display,gc,WhitePix);
+      XDrawImageString(display,pix,gc,xp,yp,Text,len);
+      strcpy(Text,"High ");
+      len    = strlen(Text);
+      width  = XTextWidth(FixFont,Text,len);
+      xp     = WinXSize/2 - AMPLI_COLS * 10 - width;
+      yp     = WinYSize - FixFontHeight + 10;
+      r      = FixFontHeight-2;
+      XSetForeground(display,gc,WhitePix);
+      if (xp > 0) XDrawImageString(display,pix,gc,xp,yp,Text,len);
+      strcpy(Text," Low");
+      len    = strlen(Text);
+      width  = XTextWidth(FixFont,Text,len);
+      xp     = WinXSize/2 + AMPLI_COLS * 10;
+      if (xp + width < WinXSize) XDrawImageString(display,pix,gc,xp,yp,Text,len);
+      yp     = InYSize;
+      XDrawLine(display,pix,gc,OffX,yp,WinXSize,yp);
+      yp     = WinYSize - r -  FixFont->descent;
+
+      for (i=0 ; i < AMPLI_COLS ; ++i)
+      {
+         xp = WinXSize/2 - AMPLI_COLS * 10 + i * 20;
+         XSetForeground(display,gc,AmpliStrip[AMPLI_COLS-i-1].pixel);
+         XFillArc(display,pix,gc,xp,yp,r,r,0,360*64);
+      }
+   }
+
    if (SizeChanged && PicType == ISOTS) /* Timeseries Caption */
    {
       XSetForeground(display,gc,BlackPix);
@@ -4593,7 +4663,7 @@ void iso(int w,int PicType,REAL *field,int dimx,int dimy,int dimz,int pal)
       XSetBackground(display,gc,BlackPix);
    }
 
-   if (PicType == ISOTS) 
+   if (PicType == ISOTS)
    {
       if (TSxp[w] == NULL) TSxp[w] = SizeAlloc(DimX * DimY , sizeof(XPoint),"TSxp");
       if (Dmin[w] == NULL) Dmin[w] = FloatAlloc(DimY ,"Dmin");
@@ -4621,7 +4691,7 @@ void iso(int w,int PicType,REAL *field,int dimx,int dimy,int dimz,int pal)
          o = Dmin[w][j];
          if ((Dmax[w][j] - Dmin[w][j]) > 1.0e-20) f = (InYSize-2) / (Dmax[w][j] - Dmin[w][j]);
          else f = 1.0;
-   
+
          for (i=1 ; i < DimX ; ++i)
          {
             TSxp[w][i].x = VGAX * i;
@@ -4631,7 +4701,7 @@ void iso(int w,int PicType,REAL *field,int dimx,int dimy,int dimz,int pal)
       }
    }
 
-   if (PicType == ISOTAB) 
+   if (PicType == ISOTAB)
    {
       XSetForeground(display,gc,BlackPix);
       XSetBackground(display,gc,BlackPix);
@@ -4783,26 +4853,26 @@ void iso(int w,int PicType,REAL *field,int dimx,int dimy,int dimz,int pal)
    {
       IsoAreas(Cstrip);
    }
-   
+
    // amplitudes of coeeficients of spherical harmonics
 
    if (PicType == ISOSH)
    {
-      SH_Amplitudes(w);     
+      SH_Amplitudes(w);
    }
 
    // amplitudes of rectangular array
 
    if (PicType == ISOAMP)
    {
-      Amplitudes(w);     
+      FC_Amplitudes(w);
    }
 
    // line plot
 
    if (PicType == LINE)
    {
-      LinePlot(w);     
+      LinePlot(w);
    }
 
    if (Grid && PicType == ISOCS) ShowGridCS();