Index: /palm/trunk/DOC/misc/plant_canopy_model.tex
===================================================================
--- /palm/trunk/DOC/misc/plant_canopy_model.tex	(revision 187)
+++ /palm/trunk/DOC/misc/plant_canopy_model.tex	(revision 187)
@@ -0,0 +1,109 @@
+\documentclass[11pt,a4paper,titlepage]{scrreprt}
+\documentstyle[Flow]
+\usepackage{graphics}
+\usepackage{german}
+\usepackage{chimuk}
+\usepackage{bibgerm}
+\usepackage{a4wide}
+\usepackage{amsmath}
+\usepackage{flafter}
+\usepackage[dvips]{epsfig}
+\usepackage{texdraw}
+\usepackage{supertabular}
+\usepackage{longtable}
+\usepackage{scrpage}
+\usepackage{fancyheadings}
+\usepackage{flow}
+\frenchspacing
+\sloppy
+\pagestyle{fancyplain}
+\addtolength{\headheight}{\baselineskip}
+\renewcommand{\chaptermark}[1]{\markboth{\thechapter~~#1}{}}
+\renewcommand{\sectionmark}[1]{\markright{\thesection\ #1}}
+\lhead{\fancyplain{}{\bfseries\leftmark}}
+\rhead{}
+%\begin{titlepage}
+%\author{Gerald Steinfeld}
+%\title{Prandtl layer parameterisation in PALM}
+%\end{titlepage}
+\begin{document}
+%\maketitle
+%\tableofcontents
+\pagebreak
+
+\chapter{Plant canopy model in PALM}
+The basic equations used in the LES are the Navier-Stokes equations using the Boussinesq approximation (conservation of momentum; equation \ref{NSG}), 
+the equation of continuity (conservation of mass; equation \ref{KON}), the first law of thermodynamics (conservation of energy; equation \ref{1HS}) and 
+the conservation equation for a passive scalar (equation \ref{CPS}). All equations are filtered (in PALM by averaging over a discrete grid volume) in 
+order to eliminate small-scale perturbations. Thus, the set of equations on that the LES model PALM is based is formed by the following equations:  
+
+\begin{equation} \label{NSG}
+\frac{\partial \overline{u}_i}{\partial t} = - \frac{\partial}{\partial x_k} \overline{u}_k \overline{u}_i - \frac{1}{\rho_0} \frac{\partial \overline{p^{\ast}}}{\partial x_i} - \left ( \epsilon_{ijk} f_j \overline{u}_k - \epsilon_{i3k} f_3 u_{g_k} \right ) + g \frac{\overline{\Theta^{\ast}}}{\Theta_0} \delta_{i3} - \frac{\partial}{\partial x_k} \overline{u_i'u_k'} - c_d a U \overline{u}_i
+\end{equation}
+
+\begin{equation} \label{KON} 
+\frac{\partial \overline{u}}{\partial x_i} = 0
+\end{equation}
+
+\begin{equation} \label{1HS}
+\frac{\partial \overline{\Theta}}{\partial t} = - \frac {\partial}{\partial x_k} \overline{u}_k \overline{\Theta} - \frac{\partial}{\partial x_k} \overline{u_k'\Theta'} + S
+\end{equation}
+
+\begin{equation} \label{CPS}
+\frac{\partial \overline s}{\partial t} = - \frac{\partial}{\partial x_k} \overline{u}_k \overline{s} - \frac{\partial}{\partial x_k} \overline{u_k's'} - c_s a U \left ( \overline{s} - s_l \right )
+\end{equation}
+
+The overline denotes filtered variables, while the prime denotes sub-gridscale, local spatial fluctuations at time t that are not resolved by the model; 
+the star denotes a departure of the basic state that is marked by a subscript $0$ and that varies only 
+with height; $u_i$ is the velocity in the $x_i$-direction; $\rho$ is the air density; $p$ is the pressure; $f_i=(0,2\Omega cos(\phi),2\Omega sin(\phi))$ 
+is the Coriolis parameter with the angular velocity of the earth $\Omega$; ${u_g}_i$ denotes the components of the geostrophic wind; $g$ is the 
+acceleration due to gravity; $\Theta$ is the potential temperature and $s$ denotes the concentration of a passive scalar.
+The sub-gridscale fluxes in equations \ref{NSG}, \ref{1HS} and \ref{CPS} are parameterised by a 1.5-order closure model using resolved-scale variables 
+and the sub-gridscale turbulent kinetic energy $\overline{e}$. Thus, an additional prognostic equation has to be solved for the sub-gridscale turbulent 
+kinetic energy $\overline{e}$:     
+
+\begin{equation} \label{TKE}
+\frac{\partial \overline e}{\partial t} = - \frac{\partial}{\partial x_k} \overline{u}_k \overline{e} - \tau_{ij} \frac{\partial \overline{u}_i}{\partial x_j} + \frac{g}{\Theta_0} \overline{u_3'\Theta_j'} - \frac{\partial}{\partial x_j} \left ( \overline{u_j' \left ( e + \frac{p'}{\rho_0} \right ) } \right ) - \epsilon - 2 c_d a U \overline{e}.
+\end{equation}
+
+Here, $\tau_{ij}$ denotes the sub-gridscale stress tensor. The final terms in equations \ref{NSG}, \ref{1HS}, \ref{CPS} and \ref{TKE} represent the sink 
+or source of momentum, heat and the passive scalar due to canopy elements, respectively. Thus, these terms do only occur within the plant canopy, while 
+they vanish outside the plant canopy (in PALM this distinction between canopy layer and the remaining area is realised by the specification of a height 
+index that contains the information up to which vertical grid level the plant canopy extends; the additional terms for the plant canopy layer are 
+only evaluated for grid levels smaller or equal that plant canopy height index). The additional, canopy related source and sink sources, are considered in 
+PALM in the same way as in the models of WATANABE (2004, Boundary-Layer Meteorol., 112, 307-341) and SHAW AND SCHUMANN 
+(1992, Boundary-Layer Meteorol., 61, 47-64).   
+$c_d$ and $c_s$ are the leaf drag coefficient and the scalar exchange 
+coefficient for a leaf, respectively; $a$ is the leaf area density; $U$ is the instantaneous local wind speed, defined as 
+\begin{equation}
+U = \left (  \overline{u}^2 + \overline{v}^2 + \overline{w}^2 \right )^{\frac{1}{2}},
+\end{equation}
+and $s_l$ is the scalar concentration at a leaf surface. 
+The heat source distribution $S(z)$ in convective conditions has been given by BROWN AND COVEY (1966, Agric. Meteorol., 3, 73-96) and SHAW AND SCHUMANN 
+(1992). The assumption that leads to the heat source distribution is that short-wave solar radiation penetrates the plant canopy and warms the foliage 
+which, in turn, warms the air in contact with it. The heat source is distributed in such a way that the heat flux $Q$ at the top of the canopy, $z=h$, is 
+prescribed and that the heat flux within the canopy follows a declining exponential function of the downward cumulative leaf area index. The decline of 
+heat flux with decreasing height is similar to the decline of net radiation that has been described by BROWN AND COVEY (1966). The heat source $S(z)$ that 
+is included in equation \ref{1HS} is the vertical derivative of the upward kinematic vertical heat flux $Q$   
+\begin{equation}
+S = \frac{\partial Q}{\partial z}.
+\end{equation}
+According to the explanations given above the upward kinematic vertical heat flux $Q$ is evaluated as
+\begin{equation} \label{kvh}
+Q(z) = Q(h) exp(- \alpha F(z)),
+\end{equation}
+where $F$ is the downward cumulative leaf area index (non-dimensional) determined by integration over the leaf area density $a$ 
+\begin{equation}
+F = \int_z^h a dz
+\end{equation}
+and $\alpha$ is an extinction coefficient that is set by default to 0.6 in PALM.
+
+The final term in equation \ref{TKE} represents the sink for sub-gridscale turbulent kinetic energy 
+due to the rapid dissipation of wake turbulence formed in the lee of canopy elements. This term is included 
+on the assumption that individual wake motions behind canopy elements are of even smaller scale than those 
+making up the bulk of sub-gridscale kinetic energy. For simulations in that the grid length is of the 
+order of the scale of wake turbulence modifications within equation \ref{TKE} might be necessary. KANDA 
+AND HINO (1994, Boundary-Layer Meteorol., 68, 237-257) made a suggestion how equation \ref{TKE} has to be 
+modified for such cases. 
+
+\end{document}
Index: /palm/trunk/DOC/misc/prandtl_layer.tex
===================================================================
--- /palm/trunk/DOC/misc/prandtl_layer.tex	(revision 187)
+++ /palm/trunk/DOC/misc/prandtl_layer.tex	(revision 187)
@@ -0,0 +1,249 @@
+\documentclass[11pt,a4paper,titlepage]{scrreprt}
+\documentstyle[Flow]
+\usepackage{graphics}
+\usepackage{german}
+\usepackage{chimuk}
+\usepackage{bibgerm}
+\usepackage{a4wide}
+\usepackage{amsmath}
+\usepackage{flafter}
+\usepackage[dvips]{epsfig}
+\usepackage{texdraw}
+\usepackage{supertabular}
+\usepackage{longtable}
+\usepackage{scrpage}
+\usepackage{fancyheadings}
+\usepackage{flow}
+\frenchspacing
+\sloppy
+\pagestyle{fancyplain}
+\addtolength{\headheight}{\baselineskip}
+\renewcommand{\chaptermark}[1]{\markboth{\thechapter~~#1}{}}
+\renewcommand{\sectionmark}[1]{\markright{\thesection\ #1}}
+\lhead{\fancyplain{}{\bfseries\leftmark}}
+\rhead{}
+%\begin{titlepage}
+%\author{Gerald Steinfeld}
+%\title{Prandtl layer parameterisation in PALM}
+%\end{titlepage}
+\begin{document}
+%\maketitle
+%\tableofcontents
+\pagebreak
+
+\chapter{Prandtl layer parameterisation in PALM}
+
+The friction velocity $u_{\ast}$ is a velocity scale that is defined by the relation 
+\begin{equation}
+u_{\ast}=\left ( \left | \tau / \overline{\rho} \right | \right )^{\frac{1}{2}}, 
+\end{equation}
+where $\tau$ is the Reynolds stress and $\rho$ is the air density.
+Using the surface kinematic momentum fluxes in the x and y directions 
+$\left ( -\overline{u'w'}, -\overline{v'w'} \right )$ to represent 
+surface stress, the friction velocity can be written as   
+\begin{equation}
+u_{\ast}=\left ( (-\overline{u'w'})^2 + (-\overline{v'w'})^2 \right )^{\frac{1}{4}}.
+\end{equation}
+Based on the definition of the friction velocity $u_{\ast}$ the vertical turbulent 
+momentum flux $\overline{u'w'}$ can be determined from 
+\begin{equation}
+{u^2_{\ast}} \cos \left ( \alpha_0 \right ) = - \overline{u'w'},
+\end{equation}
+while the vertical turbulent momentum flux $\overline{v'w'}$ can be determined 
+from  
+\begin{equation}
+{u^2_{\ast}} \sin \left ( \alpha_0 \right ) = - \overline{v'w'}.
+\end{equation}
+The angle $\alpha_0$, that is assumed to be constant in the Prandtl layer, is the 
+angle between the x-direction and the direction of the mean horizontal wind and can 
+be evaluated by   
+\begin{equation} \label{win}
+\alpha_0 = \arctan \left ( \frac{\overline{v}(z_p)}{\overline{u}(z_p)} \right ).
+\end{equation}
+According to the similarity theory of Monin and Obukhov the following relationship for 
+the profile of the mean horizontal wind is valid
+\begin{equation} \label{hmpr}
+\begin{split}
+\frac{\partial \left | \overline{\vec{v}} \right |}{\partial z} &= \frac{u_{\ast}}{\kappa z} \phi_m \left ( \frac{z}{L} \right ) \\
+                                                                &= u_{\ast} \frac{1}{\kappa z} \phi_m \left ( \frac{z}{L} \right ).
+\end{split}  
+\end{equation}
+In equation \ref{hmpr} $L$ denotes the Monin Obukhov length and $\kappa$ denotes the von Karman constant, while $\phi_m$ denotes the profile 
+or Dyer-Businger functions for momentum:
+\begin{equation}
+\phi_m = 
+\begin{cases}
+1+5 \text{Rif} & \text{if $\text{Rif} > 0$}, \\
+1 & \text{if $\text{Rif} = 0$}, \\
+\left ( 1 - 16 \text{Rif} \right )^{-\frac{1}{4}} & \text{if $\text{Rif} < 0$}. 
+\end{cases}
+\end{equation} 
+Rif denotes the dimensionless Richardson flux number.
+
+By integrating equation \ref{hmpr} over $z$ from $z_0$ to a height $z$ the following relationship for the friction velocity can be deduced:    
+\begin{equation} \label{usb}
+u_{\ast} =
+\begin{cases}
+\frac{\left | \overline{\vec{v}} \right | \kappa}{\left ( \ln \left ( \frac{z}{z_0} \right ) + 5 \text{Rif} \left ( \frac{z-z_0}{z} \right ) \right )} & \text{if $\text{Rif} \ge 0$}, \\
+%\frac{\left | \overline{\vec{v}} \right | \kappa}{\left ( ln \left ( \frac{1+B}{1-B} \frac{1-A}{1+A} \right ) + 2 \left ( arctan(B) - arctan(A) \right ) \right )} & \text{if $Rif<0$}
+\frac{\left | \overline{\vec{v}} \right | \kappa}{\ln \left ( \frac{z}{z_0} \right ) - \ln \left ( \frac{ \left ( 1+A \right )^2  \left ( 1+A^2 \right ) }{ \left ( 1+B \right )^2  \left ( 1+B^2 \right ) }\right ) + 2 \left ( \arctan(A) - \arctan(B) \right )} & \text{if $\text{Rif}<0$}
+\end{cases} 
+\end{equation}
+with
+\begin{equation}
+%A=\left ( 1 - 16 Rif \right )^{-\frac{1}{4}}
+A=\left ( 1 - 16 \text{Rif} \right )^{\frac{1}{4}} 
+\end{equation}
+and
+\begin{equation}
+%B=\left ( 1 - 16 Rif \frac{z_0}{z} \right )^{-\frac{1}{4}}.
+B=\left ( 1 - 16 \text{Rif} \frac{z_0}{z} \right )^{\frac{1}{4}}.
+\end{equation}
+In fact, equation \ref{usb} is used in PALM to determine the friction velocity $u_{\ast}$.
+
+The following paragraph is a short digression dealing with the integration of equation \ref{hmpr} that finally leads to equation \ref{usb}. The integration 
+is only shown for the profile function for unstable stratification. According to PAULSON (1970) the general result of an integration of $\phi=\frac{\kappa z'}{u_{\ast}} 
+\frac{\partial \overline{u}}{\partial z'}$ over $z'$ from $z_0$ to $z$ is
+\begin{equation}
+\overline{u}(z) = \frac{u_{\ast}}{\kappa} \left [ \ln \left ( \frac{z}{z_0} \right ) - \Psi \right ] 
+\end{equation}
+with
+\begin{equation}
+\Psi = \int_{\frac{z_0}{L}}^{\frac{z}{L}}d \left ( \frac{z'}{L} \right ) \frac{1-\phi \left ( \frac{z'}{L} \right )}{\frac{z'}{L}}.
+\end{equation}  
+Applying the profile function for unstable stratification, $\phi=\left ( 1 - 16 \frac{z}{L} \right )^{-\frac{1}{4}}$, $\Psi$ can be determined as follows:
+\begin{equation}
+\begin{split}
+\Psi &= \int_{\frac{z_0}{L}}^{\frac{z}{L}}d \left ( \frac{z'}{L} \right ) \frac{1- \left ( 1 - 16 \frac{z'}{L} \right )^{-\frac{1}{4}} }{\frac{z'}{L}} \text{ | substitution: } x=\frac{1}{\phi \left ( \frac{z'}{L} \right ) } \\
+     &= \int_{\left ( 1-16\frac{z_0}{L} \right )^{\frac{1}{4}}}^{\left ( 1-16\frac{z}{L} \right )^{\frac{1}{4}}} dx \left ( 4 \frac{x^2 - x^3}{1 - x^4} \right ) \\
+     &= \int_{\left ( 1-16\frac{z_0}{L} \right )^{\frac{1}{4}}}^{\left ( 1-16\frac{z}{L} \right )^{\frac{1}{4}}} dx \left ( 2 \left ( \frac{1}{1+x} + \frac{x}{1+x^2} - \frac{1}{1+x^2} \right ) \right ) \\
+     &= \left [ 2 \ln \left ( \frac{1+x}{2} \right ) + \ln \left ( \frac{1+x^2}{2} \right ) -2 \arctan(x) \right ]_{\left ( 1-16\frac{z_0}{L} \right )^{\frac{1}{4}}=B}^{\left ( 1-16\frac{z}{L} \right )^{\frac{1}{4}}=A} \\
+     &= \ln \left ( \frac{\left ( 1 + A \right )^2}{\left ( 1 + B \right )^2} \frac{1+A^2}{1+B^2} \right ) - 2 \left ( \arctan(A) - \arctan(B) \right )
+\end{split}
+\end{equation} 
+
+Making use of equation \ref{win} and \ref{hmpr}, it is also possible to deduce a relationship for the profile of the mean u-component of the wind velocity
+\begin{equation} \label{upr}
+\begin{split}
+\frac{\partial \overline{u} }{\partial z} &= \frac{\partial \left | \overline{\vec{v}} \right |}{\partial z} \cos \left ( \alpha_0 \right ) \\
+                                          &= \frac{u_{\ast}}{\kappa z} \phi_m \left ( \frac{z}{L} \right ) \cos \left ( \alpha_0 \right ) \\
+                                          &= \frac{1}{u_{\ast}} {u^2_{\ast}} \cos \left ( \alpha_0 \right ) \frac{1}{\kappa z} \phi_m \left ( \frac{z}{L} \right ) \\
+                                          &= \frac{-\overline{u'w'}}{u_{\ast}} \frac{1}{\kappa z} \phi_m \left ( \frac{z}{L} \right )
+\end{split} 
+\end{equation}
+and accordingly a relationship for the profile of the mean v-component of the wind velocity
+\begin{equation} \label{vpr}
+\begin{split}
+\frac{\partial \overline{v} }{\partial z} &= \frac{\partial \left | \overline{\vec{v}} \right |}{\partial z} \sin \left ( \alpha_0 \right ) \\
+                                          &= \frac{u_{\ast}}{\kappa z} \phi_m \left ( \frac{z}{L} \right ) \sin \left ( \alpha_0 \right ) \\
+                                          &= \frac{1}{u_{\ast}} {u^2_{\ast}} \sin \left ( \alpha_0 \right ) \frac{1}{\kappa z} \phi_m \left ( \frac{z}{L} \right ) \\
+                                          &= \frac{-\overline{v'w'}}{u_{\ast}} \frac{1}{\kappa z} \phi_m \left ( \frac{z}{L} \right ).
+\end{split} 
+\end{equation}
+As the right-hand sides of equation \ref{upr} and \ref{vpr} differ only by prefactors that are (assumed to be) constant within the Prandtl layer from the 
+right-hand side of equation \ref{hmpr}, we can directly make use of the integration that led to equation \ref{usb} in order to obtain   
+\begin{equation} \label{uswsb}
+C=\frac{-\overline{u'w'}}{u_{\ast}} =
+\begin{cases}
+\frac{\overline{u} \kappa}{\left ( \ln \left ( \frac{z}{z_0} \right ) + 5 \text{Rif} \left ( \frac{z-z_0}{z} \right ) \right )} & \text{if $\text{Rif} \ge 0$}, \\
+%\frac{\overline{u} \kappa}{\left ( ln \left ( \frac{1+B}{1-B} \frac{1-A}{1+A} \right ) + 2 \left ( arctan(B) - arctan(A) \right ) \right )} & \text{if $Rif<0$}.
+\frac{ \overline{u} \kappa}{\ln \left ( \frac{z}{z_0} \right ) - \ln \left ( \frac{ \left ( 1+A \right )^2  \left ( 1+A^2 \right ) }{ \left ( 1+B \right )^2  \left ( 1+B^2 \right ) }\right ) + 2 \left ( \arctan(A) - \arctan(B) \right )} & \text{if $\text{Rif}<0$}.
+
+\end{cases} 
+\end{equation}
+and
+\begin{equation} \label{vswsb}
+D=\frac{-\overline{v'w'}}{u_{\ast}} =
+\begin{cases}
+\frac{\overline{v} \kappa}{\left ( \ln \left ( \frac{z}{z_0} \right ) + 5 \text{Rif} \left ( \frac{z-z_0}{z} \right ) \right )} & \text{if $\text{Rif} \ge 0$}, \\
+%\frac{\overline{v} \kappa}{\left ( ln \left ( \frac{1+B}{1-B} \frac{1-A}{1+A} \right ) + 2 \left ( arctan(B) - arctan(A) \right ) \right )} & \text{if $Rif<0$}.
+\frac{\overline{v} \kappa}{\ln \left ( \frac{z}{z_0} \right ) - \ln \left ( \frac{ \left ( 1+A \right )^2  \left ( 1+A^2 \right ) }{ \left ( 1+B \right )^2  \left ( 1+B^2 \right ) }\right ) + 2 \left ( \arctan(A) - \arctan(B) \right )} & \text{if $\text{Rif}<0$}.
+
+\end{cases} 
+\end{equation}
+Both equations, \ref{uswsb} and \ref{vswsb}, are used in PALM. In order to get an information on the turbulent vertical momentum fluxes within the Prandtl layer, 
+equation \ref{uswsb} and equation \ref{vswsb} need only to be multiplied by $-1$ and the friction velocity $u_{\ast}$, so that finally the followings two steps 
+have to be executed in PALM: 
+\begin{equation}
+\overline{u'w'} = -C u_{\ast}
+\end{equation}
+and
+\begin{equation}
+\overline{v'w'} = -D u_{\ast}.
+\end{equation}
+
+In case that no near-surface heat flux has been prescribed by the user of PALM, the near-surface heat flux $\overline{w'\Theta'}_0$ is evaluated from
+\begin{equation} \label{hfe}
+\overline{w'\Theta'} = - u_{\ast} \theta_{\ast}.
+\end{equation}
+Here, $\theta_{\ast}$ is the so-called characteristic temperature for the Prandtl layer. In case of no preset near-surface heat flux the characteristic 
+temperature $\theta_{\ast}$ is determined in PALM from the integrated version of the profile function for the potential temperature. According to the 
+similarity theory of Monin and Obukhov the following relationship for the vertical gradient of the potential temperature is valid:
+\begin{equation} \label{hmpt}
+\frac{\partial \overline{\Theta}}{\partial z} = \frac{\theta_{\ast}}{\kappa z} \phi_h.
+\end{equation}  
+In equation \ref{hmpt} $\phi_h$ denotes the profile or Dyer-Businger functions for temperature:
+\begin{equation}
+\phi_h = 
+\begin{cases}
+1+5 \text{Rif} & \text{if $\text{Rif} > 0$}, \\
+1 & \text{if $\text{Rif} = 0$}, \\
+\left ( 1 - 16 \text{Rif} \right )^{-\frac{1}{2}} & \text{if $\text{Rif} < 0$}. 
+\end{cases}
+\end{equation} 
+By integrating equation \ref{hmpt} over $z$ from $z_0$ to a height $z$ the following relationship for the characteristic temperature in the Prandtl 
+layer can be deduced:    
+\begin{equation} \label{tsb}
+\theta_{\ast} =
+\begin{cases}
+\frac{\kappa \left ( \overline{\Theta}(z) - \overline{\Theta}(z_0) \right )}{\ln \left ( \frac{z}{z_0} \right ) + 5 \text{Rif} \left ( \frac{z-z_0}{z} \right )} & \text{if $\text{Rif} \ge 0$}, \\
+\frac{\kappa \left ( \overline{\Theta}(z) - \overline{\Theta}(z_0) \right )}{\ln \left ( \frac{z}{z_0} - 2 \ln \left ( \frac{1+A}{1+B} \right )\right )} & \text{if $\text{Rif}<0$}
+\end{cases} 
+\end{equation}
+with
+\begin{equation}
+%A=\left ( 1 - 16 \text{Rif} \right )^{-\frac{1}{4}}
+A=\sqrt{1 - 16 \text{Rif}} 
+\end{equation}
+and
+\begin{equation}
+%B=\left ( 1 - 16 \text{Rif} \frac{z_0}{z} \right )^{-\frac{1}{4}}.
+B=\sqrt{1 - 16 \text{Rif} \frac{z_0}{z}}.
+\end{equation}
+Note, that the temperature at the height of the roughness length that is required for the evaluation of $\theta_{\ast}$ in equation \ref{tsb} is saved 
+in PALM at the first vertical grid level of the temperature with index k=0. The height that is assigned to this vertical grid level is -0.5$\Delta z$, 
+where $\Delta z$ is the grid length in the vertical direction.   
+In case of a near-surface heat flux $\overline{w'\Theta'}_0$ that has been prescribed by the user, the evaluation of $\theta_{\ast}$ is not needed for 
+the integration of the model. Instead $\theta_{\ast}$ is only determined for the evaluation of statistics of the turbulent flow. In that case it is 
+simply derived by a transformation of equation \ref{hfe} into  
+\begin{equation}
+\theta_{\ast} = - \frac{\overline{w'\Theta'}_0}{u_{\ast}}.
+\end{equation}    
+In case that PALM is run in its moist version, the evaluation of the characteristic humidity for the Prandtl layer is evaluated in a way correspondent to 
+that of the determination of the characteristic temperature.
+
+
+  
+  
+ 
+        
+
+ 
+
+ 
+
+
+ 
+
+
+
+
+
+
+
+ 
+  
+
+ 
+
+\end{document}
Index: /palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 186)
+++ /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 187)
@@ -42,4 +42,8 @@
 Changed:
 -------
+Modification of the integrated version of the profile function for momentum 
+for unstable stratification; more consistent flux definitions. (wall_fluxes, 
+production_e)
+
 Strict grid matching along z is not needed for mg-solver. (check_parameters)
 
@@ -79,9 +83,14 @@
 (init_1d_model)
 
-advec_s_ups, advec_u_ups, advec_v_ups, advec_w_ups, calc_spectra, check_open, check_parameters, cpu_statistics, init_1d_model, init_3d_model, modules, palm, parin, poisfft, read_var_list, read_3d_binary, transpose, write_var_list, write_3d_binary
+advec_s_ups, advec_u_ups, advec_v_ups, advec_w_ups, calc_spectra, check_open, check_parameters, cpu_statistics, init_1d_model, init_3d_model, modules, palm, parin, poisfft, production_e, read_var_list, read_3d_binary, transpose, wall_fluxes, write_var_list, write_3d_binary
 
 
 Errors:
 ------
+Bugfix: change definition of us_wall from 1D to 2D. Tests showed that this decreases 
+u* by some 10% and increases TKE and momentum fluxes by some 10% because friction 
+was underestimated and momentum fluxes were wrongly calculated due to the bug. 
+(prandtl_fluxes, wall_fluxes)
+
 Bugfix: calculation of horizontal fluxes at vertical walls (diffusion_s)
 
@@ -106,20 +115,27 @@
 
 Introduce prefix_chr to ensure unique dvrp_file path.
+
 small bugfixes for user_interface sample code (comments):
 - initialize ustvst with 0.0 as it is now computed only until nxr and nyn
 - two ALLOCATE statements moved from user_read_restart_data back to user_init
 - remove 'READ (13) u2_av' statement in user_read_restart_data
+
 Bugfix: remove IF statement in plant_canopy_model_ij (plant_canopy_model)
+
 Bugfix: divide sums(k,8) (e) and sums(k,34) (e*) by ngp_2dh_s_inner(k,sr) 
 (like other scalars) (flow_statistics)
+
 Bugfix: dopr_time_count was written on the binary file, which caused that
 NetCDF files newly created by restart files (no append of existing files!)
 contained uneccessary time levels. (read_3d_binary, write_3d_binary)
+
 Bugfix: extra '*' removed in user_statistics sample code (user_interface)
+
 Bugfix: a stop command was missing in some cases of the parallel branch (local_stop)
+
 Bugfix in volume flow control for non-cyclic boundary conditions (pres)
+
 Bugfix: misplaced #endif directives (combine_plot_fields)
 
+check_parameters, diffusion_s, flow_statistics, init_dvrp, init_3d_model, local_stop, plant_canopy_model, poismg, prandtl_fluxes, pres, read_3d_binary, user_interface, wall_fluxes, write_3d_binary
 
-check_parameters, diffusion_s, flow_statistics, init_dvrp, init_3d_model, local_stop, plant_canopy_model, poismg, pres, read_3d_binary, user_interface, write_3d_binary
-
Index: /palm/trunk/SOURCE/prandtl_fluxes.f90
===================================================================
--- /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 186)
+++ /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 187)
@@ -4,5 +4,20 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Bugfix: change definition of us_wall from 1D to 2D:
+!         Modification of the evaluation of the vertical turbulent momentum 
+!         fluxes u'w' and v'w'; the first usws that is computed corresponds 
+!         to -u'w'/u* and not as priorily assumed to (-u'w')**0.5, the first 
+!         vsws that is computed corresponds to -v'w'/u* and not as priorily
+!         assumed to (-v'w')**0.5. Therefore, the intermediate result for 
+!         usws has to be multiplied by -u* instead by itself in order to 
+!         get u'w'. Accordingly, the intermediate result for vsws has to be 
+!         multiplied by -u* instead by itself in order to get v'w'. As u* 
+!         is calculated for the position of a scalar an additional 
+!         interpolation of u* to the position of u and v, respectively, 
+!         is necessary. As u* is not determined for the ghost points on 
+!         each PE, an additional exchange of information from neighbouring 
+!         PEs is necessary.!  
+! Change: Modification of the integrated version of the profile function for 
+!         momentum for unstable stratification
 !
 ! Former revisions:
@@ -81,19 +96,9 @@
 !--             Unstable stratification
                 a = SQRT( 1.0 - 16.0 * rif(j,i) )
-                b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
-!
-!--             If a borderline case occurs, the formula for stable
-!--             stratification must be used anyway, or else a zero division
-!--             would occur in the argument of the logarithm
-                IF ( a == 1.0  .OR.  b == 1.0 )  THEN
-                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
-                                     LOG( z_p / z0(j,i) ) +                   &
-                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
-                                                                   )
-                ELSE
-                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
-                                 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
-                                                                   )
-                ENDIF
+                b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
+
+                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (         &
+                          LOG( z_p / z0(j,i) ) -                           &
+                          2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) )
              ENDIF
 
@@ -168,25 +173,21 @@
 !
 !--          Unstable stratification
-             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
-             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
-!
-!--          If a borderline case occurs, the formula for stable stratification
-!--          must be used anyway, or else a zero division would occur in the 
-!--          argument of the logarithm.
-             IF ( a == 1.0  .OR.  b == 1.0 )  THEN
-                us(j,i) = kappa * uv_total / (                                &
-                                     LOG( z_p / z0(j,i) ) +                   &
-                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
-                                             )
-             ELSE
-                us(j,i) = kappa * uv_total / (                               &
-                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
-                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
-                                             )
-             ENDIF
+             a = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
+             b = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
+
+             us(j,i) = kappa * uv_total / (                                  &
+                       LOG( z_p / z0(j,i) ) -                                &
+                       LOG( ( 1.0 + a )**2 * ( 1.0 + a**2 ) / (              &
+                            ( 1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +          &
+                            2.0 * ( ATAN( a ) - ATAN( b ) )                  &
+                                           )
           ENDIF
        ENDDO
     ENDDO
 
+!
+!-- Values of us at ghost point locations are needed for the evaluation of usws 
+!-- and vsws.
+    CALL exchange_horiz_2d( us )
 !
 !-- Compute u'w' for the total model domain.
@@ -212,23 +213,15 @@
 !
 !--          Unstable stratification
-             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
-             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
-!
-!--          If a borderline case occurs, the formula for stable stratification
-!--          must be used anyway, or else a zero division would occur in the 
-!--          argument of the logarithm. 
-             IF ( a == 1.0  .OR.  B == 1.0 )  THEN
-                usws(j,i) = kappa * u(k+1,j,i) / (                           &
-                                        LOG( z_p / z0(j,i) ) +               &
-                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
+             a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
+             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
+
+             usws(j,i) = kappa * u(k+1,j,i) / (                           &
+                         LOG( z_p / z0(j,i) ) -                           &
+                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
+                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
+                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
                                                  )
-             ELSE
-                usws(j,i) = kappa * u(k+1,j,i) / (                           &
-                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
-                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
-                                                 )
-             ENDIF
           ENDIF
-          usws(j,i) = -usws(j,i) * ABS( usws(j,i) )
+          usws(j,i) = -usws(j,i) * 0.5 * ( us(j,i-1) + us(j,i) )
        ENDDO
     ENDDO
@@ -257,23 +250,15 @@
 !
 !--          Unstable stratification
-             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
-             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
-!
-!--          If a borderline case occurs, the formula for stable stratification
-!--          must be used anyway, or else a zero division would occur in the 
-!--          argument of the logarithm.
-             IF ( a == 1.0  .OR.  b == 1.0 )  THEN
-                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
-                                        LOG( z_p / z0(j,i) ) +               &
-                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
+             a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
+             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
+
+             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
+                         LOG( z_p / z0(j,i) ) -                           &
+                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
+                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
+                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
                                                  )
-             ELSE
-                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
-                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
-                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
-                                                 )
-             ENDIF
           ENDIF
-          vsws(j,i) = -vsws(j,i) * ABS( vsws(j,i) )
+          vsws(j,i) = -vsws(j,i) * 0.5 * ( us(j-1,i) + us(j,i) )
        ENDDO
     ENDDO
@@ -317,20 +302,10 @@
 !
 !--                Unstable stratification
-                   a = SQRT( 1.0 - 16.0 * rif(j,i) )
-                   b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
-!
-!--                If a borderline case occurs, the formula for stable
-!--                stratification must be used anyway, or else a zero division
-!--                would occur in the argument of the logarithm. 
-                   IF ( a == 1.0  .OR.  b == 1.0 )  THEN
-                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
-                                     LOG( z_p / z0(j,i) ) +                   &
-                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
-                                                                    )
-                   ELSE
-                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
-                                  LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
-                                                                    )
-                   ENDIF
+                   a = SQRT( 1.0 - 16.0 * rif(j,i) ) 
+                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p ) 
+ 
+                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (   &
+                             LOG( z_p / z0(j,i) ) -                    &
+                              2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
                 ENDIF
 
@@ -341,7 +316,6 @@
 
 !
-!-- Exchange the boundaries for u* and the momentum fluxes (fluxes only for
-!-- completeness's sake).
-    CALL exchange_horiz_2d( us )
+!-- Exchange the boundaries for the momentum fluxes (only for sake of 
+!-- completeness)
     CALL exchange_horiz_2d( usws )
     CALL exchange_horiz_2d( vsws )
Index: /palm/trunk/SOURCE/production_e.f90
===================================================================
--- /palm/trunk/SOURCE/production_e.f90	(revision 186)
+++ /palm/trunk/SOURCE/production_e.f90	(revision 187)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Change: add 'minus' sign to fluxes obtained from subroutine wall_fluxes_e for
+!         consistency with subroutine wall_fluxes
 !
 ! Former revisions:
@@ -164,10 +165,10 @@
                       CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                           usvs, 1.0, 0.0, 0.0, 0.0 )
-                      dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
-!                      dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
+                      dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
+!                      dudy = - wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
                       CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                           wsvs, 0.0, 0.0, 1.0, 0.0 )
-                      dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
-!                      dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
+                      dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
+!                      dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
                    ELSE
                       dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
@@ -180,10 +181,10 @@
                       CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                           vsus, 0.0, 1.0, 0.0, 0.0 )
-                      dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
-!                      dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
+                      dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
+!                      dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
                       CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                           wsus, 0.0, 0.0, 0.0, 1.0 )
-                      dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
-!                      dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
+                      dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
+!                      dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
                    ELSE
                       dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
@@ -220,8 +221,8 @@
 
                       IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                         dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
-!                         dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
-                         dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
-!                         dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
+                         dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
+!                         dudy = - wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
+                         dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
+!                         dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
                       ELSE
                          dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
@@ -232,8 +233,8 @@
 
                       IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                         dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
-!                         dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
-                         dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
-!                         dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
+                         dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
+!                         dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
+                         dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
+!                         dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
                       ELSE
                          dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
@@ -629,8 +630,8 @@
                 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                     usvs, 1.0, 0.0, 0.0, 0.0 )
-                dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
+                dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
                 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                     wsvs, 0.0, 0.0, 1.0, 0.0 )
-                dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
+                dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
              ELSE
                 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
@@ -643,8 +644,8 @@
                 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                     vsus, 0.0, 1.0, 0.0, 0.0 )
-                dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
+                dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
                 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                     wsus, 0.0, 0.0, 0.0, 1.0 )
-                dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
+                dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
              ELSE
                 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
@@ -679,6 +680,6 @@
 
                 IF ( wall_e_y(j,i) /= 0.0 )  THEN
-                   dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
-                   dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
+                   dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
+                   dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
                 ELSE
                    dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
@@ -689,6 +690,6 @@
 
                 IF ( wall_e_x(j,i) /= 0.0 )  THEN
-                   dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
-                   dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
+                   dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
+                   dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
                 ELSE
                    dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
Index: /palm/trunk/SOURCE/wall_fluxes.f90
===================================================================
--- /palm/trunk/SOURCE/wall_fluxes.f90	(revision 186)
+++ /palm/trunk/SOURCE/wall_fluxes.f90	(revision 187)
@@ -3,5 +3,20 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Bugfix: change definition of us_wall from 1D to 2D:
+!         Modification of the evaluation of the vertical turbulent momentum 
+!         fluxes u'w' and v'w'; the first usws that is computed corresponds 
+!         to -u'w'/u* and not as priorily assumed to (-u'w')**0.5, the first 
+!         vsws that is computed corresponds to -v'w'/u* and not as priorily
+!         assumed to (-v'w')**0.5. Therefore, the intermediate result for 
+!         usws has to be multiplied by -u* instead by itself in order to 
+!         get u'w'. Accordingly, the intermediate result for vsws has to be 
+!         multiplied by -u* instead by itself in order to get v'w'.
+!         This requires the calculation of us_wall (and vel_total, u_i, v_i, ws)
+!         also in wall_fluxes_e.
+! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed
+! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for 
+!         consistency with subroutine wall_fluxes
+! Change: Modification of the integrated version of the profile function for 
+!         momentum for unstable stratification
 !
 ! Former revisions:
@@ -69,5 +84,5 @@
 !
 !--             All subsequent variables are computed for the respective
-!--             location where the relevant variable is defined
+!--             location where the respective flux is defined.
                 DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
 
@@ -109,22 +124,13 @@
 !
 !--                   Unstable stratification
-                      h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
-                      h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ))
-
-!
-!--                   If a borderline case occurs, the formula for stable
-!--                   stratification must be used anyway, or else a zero
-!--                   division would occur in the argument of the logarithm.
-                      IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
-                         us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
-                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
-                                                       )
-                      ELSE
-                         us_wall = kappa * vel_total / (                       &
-                            LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
-                            2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
-                                                       )
-                      ENDIF
-
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      us_wall = kappa * vel_total / (                          &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                    )
                    ENDIF
 
@@ -160,27 +166,16 @@
 !
 !--                   Unstable stratification
-                      h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
-                      h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ))
-
-!
-!--                   If a borderline case occurs, the formula for stable
-!--                   stratification must be used anyway, or else a zero
-!--                   division would occur in the argument of the logarithm. 
-                      IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
-                         wall_flux(k,j,i) = kappa *                            &
-                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
-                              ( LOG( zp / z0(j,i) ) +                          &
-                                5.0 * rifs * ( zp - z0(j,i) ) / zp             &
-                              )
-                      ELSE
-                         wall_flux(k,j,i) = kappa *                            &
-                             ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /  &
-                             ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )&
-                               + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )             &
-                             )
-                      ENDIF
-
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      wall_flux(k,j,i) = kappa *                               &
+                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                                            )
                    ENDIF
-                   wall_flux(k,j,i) = -wall_flux(k,j,i) * ABS(wall_flux(k,j,i))
+                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
 
 !
@@ -226,5 +221,5 @@
 !
 !--    All subsequent variables are computed for the respective location where 
-!--    the relevant variable is defined
+!--    the respective flux is defined.
        DO  k = nzb_w, nzt_w
 
@@ -266,22 +261,13 @@
 !
 !--          Unstable stratification
-             h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
-             h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
-
-!
-!--          If a borderline case occurs, the formula for stable stratification
-!--          must be used anyway, or else a zero division would occur in the 
-!--          argument of the logarithm.
-             IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
-                us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +          &
-                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
-                                              )
-             ELSE
-                us_wall = kappa * vel_total / (                                &
-                            LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
-                            2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
-                                              )
-             ENDIF
-
+             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+             us_wall = kappa * vel_total / (                          &
+                  LOG( zp / z0(j,i) ) -                               &
+                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                           )
           ENDIF
 
@@ -315,27 +301,16 @@
 !
 !--          Unstable stratification
-             h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
-             h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
-
-!
-!--          If a borderline case occurs, the formula for stable stratification
-!--          must be used anyway, or else a zero division would occur in the 
-!--          argument of the logarithm. 
-             IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
-                wall_flux(k) = kappa *                                         &
-                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
-                              ( LOG( zp / z0(j,i) ) +                          &
-                                5.0 * rifs * ( zp - z0(j,i) ) / zp             &
-                              )
-             ELSE
-                wall_flux(k) = kappa *                                         &
-                             ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /  &
-                             ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )&
-                               + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )             &
-                             )
-             ENDIF
-
+             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+             wall_flux(k) = kappa *                               &
+                  ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
+                  LOG( zp / z0(j,i) ) -                               &
+                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                                   )
           ENDIF
-          wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
+          wall_flux(k) = -wall_flux(k) * us_wall
 
 !
@@ -371,5 +346,6 @@
 
        INTEGER ::  i, j, k, kk, wall_index
-       REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
+       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
+                   ws, zp
 
        REAL ::  rifs
@@ -388,10 +364,8 @@
              IF ( wall(j,i) /= 0.0 )  THEN
 !
-!--             All subsequent variables are computed for the respective
-!--             location where the relevant variable is defined
+!--             All subsequent variables are computed for scalar locations.
                 DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
-
-!
-!--                (1) Compute rifs
+!
+!--                (1) Compute rifs, u_i, v_i, and ws
                    IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
                       kk = nzb_diff_s_inner(j,i)-1
@@ -404,14 +378,41 @@
                                  )
 
-!
-!--                Skip (2) to (4) of wall_fluxes, because here rifs is
-!--                already available from (1)
-
+                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
+                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
+                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
+!
+!--                (2) Compute wall-parallel absolute velocity vel_total and
+!--                interpolate appropriate velocity component vel_zp.
+                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
+                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
+!
+!--                (3) Compute wall friction velocity us_wall
+                   IF ( rifs >= 0.0 )  THEN
+
+!
+!--                   Stable stratification (and neutral)
+                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
+                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                                                    )
+                   ELSE
+
+!
+!--                   Unstable stratification
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      us_wall = kappa * vel_total / (                          &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                    )
+                   ENDIF
+
+!
+!--                Skip step (4) of wall_fluxes, because here rifs is already
+!--                available from (1)
 !
 !--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
-                   vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
-                                          b * ( v(k,j,i) + v(k,j+1,i) ) +  &
-                                    (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
-                                  )
 
                    IF ( rifs >= 0.0 )  THEN
@@ -425,29 +426,15 @@
 !
 !--                   Unstable stratification
-                      h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
-                      h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ))
-
-!
-!--                   If a borderline case occurs, the formula for stable
-!--                   stratification must be used anyway, or else a zero
-!--                   division would occur in the argument of the logarithm. 
-                      IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
-                         wall_flux(k,j,i) = kappa * vel_zp /                 &
-                                        ( LOG( zp / z0(j,i) ) +              &
-                                          5.0 * rifs * ( zp - z0(j,i) ) / zp &
-                                        )
-                      ELSE
-                         wall_flux(k,j,i) = kappa * vel_zp /                   &
-                            ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
-                              + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
-                            )
-                      ENDIF
-
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                          )
                    ENDIF
-                   wall_flux(k,j,i) = wall_flux(k,j,i) * ABS( wall_flux(k,j,i) )
-
-!
-!--                Store rifs for next time step
-                   rif_wall(k,j,i,wall_index) = rifs
+                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
 
                 ENDDO
@@ -476,5 +463,6 @@
 
        INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
-       REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
+       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
+                   ws, zp
 
        REAL ::  rifs
@@ -488,10 +476,9 @@
 
 !
-!--    All subsequent variables are computed for the respective location where 
-!--    the relevant variable is defined
+!--    All subsequent variables are computed for scalar locations.
        DO  k = nzb_w, nzt_w
 
 !
-!--       (1) Compute rifs
+!--       (1) Compute rifs, u_i, v_i, and ws
           IF ( k == nzb_w )  THEN
              kk = nzb_w
@@ -504,10 +491,44 @@
                         )
 
-!
-!--       Skip (2) to (4) of wall_fluxes, because here rifs is already available
-!--       from (1)
-
+          u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
+          v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
+          ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
+!
+!--       (2) Compute wall-parallel absolute velocity vel_total and
+!--       interpolate appropriate velocity component vel_zp.
+          vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
+          vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
+!
+!--       (3) Compute wall friction velocity us_wall
+          IF ( rifs >= 0.0 )  THEN
+
+!
+!--          Stable stratification (and neutral)
+             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
+                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                                           )
+          ELSE
+
+!
+!--          Unstable stratification
+             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+             us_wall = kappa * vel_total / (                          &
+                  LOG( zp / z0(j,i) ) -                               &
+                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                           )
+          ENDIF
+
+!
+!--       Skip step (4) of wall_fluxes, because here rifs is already
+!--       available from (1)
 !
 !--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
+!--       First interpolate the velocity (this is different from 
+!--       subroutine wall_fluxes because fluxes in subroutine 
+!--       wall_fluxes_e are defined at scalar locations).
           vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
                                  b * ( v(k,j,i) + v(k,j+1,i) ) +  &
@@ -525,29 +546,15 @@
 !
 !--          Unstable stratification
-             h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
-             h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
-
-!
-!--          If a borderline case occurs, the formula for stable stratification
-!--          must be used anyway, or else a zero division would occur in the 
-!--          argument of the logarithm. 
-             IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
-                wall_flux(k) = kappa * vel_zp /                     &
-                               ( LOG( zp / z0(j,i) ) +              &
-                                 5.0 * rifs * ( zp - z0(j,i) ) / zp &
-                               )
-             ELSE
-                wall_flux(k) = kappa * vel_zp /                                &
-                            ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) &
-                              + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              &
-                            )
-             ENDIF
-
+             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+             wall_flux(k) = kappa * vel_zp / (                        &
+                  LOG( zp / z0(j,i) ) -                               &
+                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                 )
           ENDIF
-          wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) )
-
-!
-!--       Store rifs for next time step
-          rif_wall(k,j,i,wall_index) = rifs
+          wall_flux(k) = - wall_flux(k) * us_wall
 
        ENDDO
