<p><b>gill</b> 2008-04-30 11:47:21 -0600 (Wed, 30 Apr 2008)</p><p>3.0 version for Chap 2, 3, 4<br>
from Bill Skamarock<br>
<br>
M filter.tex<br>
M equation.tex<br>
M discretization.tex<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/wrf/technote/discretization.tex
===================================================================
--- trunk/wrf/technote/discretization.tex 2008-04-26 00:42:07 UTC (rev 61)
+++ trunk/wrf/technote/discretization.tex 2008-04-30 17:47:21 UTC (rev 62)
@@ -4,25 +4,28 @@
\section{Temporal Discretization}
The ARW solver uses a time-split integration scheme. Generally
-speaking, slow or low-frequency (meteorologically significant) modes
-are integrated using a third-order Runge-Kutta (RK3) time integration
+speaking, slow or low-frequency (meteorologically significant) modes are
+integrated using a third-order Runge-Kutta (RK3) time integration
scheme, while the high-frequency acoustic modes are integrated over
smaller time steps to maintain numerical stability. The horizontally
propagating acoustic modes (including the external mode present in the
mass-coordinate equations using a constant-pressure upper boundary
-condition) are integrated using a forward-backward time integration
-scheme, and vertically propagating acoustic modes and buoyancy
-oscillations are integrated using a vertically implicit scheme (using
-the acoustic time step). The time-split integration is similar to that
-first developed by \citet{klemp78} and analyzed by
-\citet{skamarock92}. The time-split RK3 scheme is described in general
-terms in \citet{wicker02}. The primary differences between the
-descriptions found in the references and the ARW implementation
-are associated with our use of the mass vertical coordinate and a
-flux-form set of equations, along with our use of perturbation
-variables for the acoustic component of the time-split integration.
-The acoustic-mode integration is cast in the form of a correction to
-the RK3 integration.
+condition) and gravity waves are integrated using a forward-backward
+time integration scheme, and vertically propagating acoustic modes and
+buoyancy oscillations are integrated using a vertically implicit scheme
+(using the acoustic time step). The time-split integration for the
+flux-form equations is described and analyzed in
+\citet{klemp_et_al_2007}. The time-splitting is similar to that first
+developed by \citet{klemp78} for leapfrog time integration and analyzed
+by \citet{skamarock92}. This time-split approach was extended to the
+RK3 scheme as described in \citet{wicker02}. The primary differences
+between the earlier implementations described in the references and the
+ARW implementation are associated with our use of the mass vertical
+coordinate and a flux-form set of equations, as described in
+\citet{klemp_et_al_2007}, along with our use of
+perturbation variables for the acoustic component of the time-split
+integration. The acoustic-mode integration is cast in the form of a
+correction to the RK3 integration.
\subsection{Runge-Kutta Time Integration Scheme}
\label{rk3_scheme}
@@ -48,7 +51,7 @@
%
</font>
<font color="gray">oindent
where $\Delta t$ is the time step for the low-frequency modes
-({\it the} model time step).
+(the model time step).
In \eqref{rk3a} -- \eqref{rk3c}, superscripts denote time levels.
This scheme is not a true Runge-Kutta scheme {\it per se} because, while it is
third-order accurate for linear equations, it is only second-order accurate
@@ -64,7 +67,7 @@
The high-frequency but meteorologically insignificant acoustic modes
would severely limit the RK3 time step $\Delta t$ in \eqref{rk3a} --
\eqref{rk3c}. To circumvent this time step limitation we use the
-approach described in \citet{wicker02}. Additionally, to increase the
+time-split approach described in \citet{wicker02}. Additionally, to increase the
accuracy of the splitting, we integrate
a perturbation form of the governing equations using smaller acoustic
time steps within the RK3 large-time-step sequence.
@@ -147,27 +150,29 @@
\label{v-small-step}
\\
\delta_\tau \mu_d''
- + m^2[\partial_x U'' + \partial_y V'']^{\tau+\Delta \tau}
-+ m \partial_\eta {\Omega''^{\tau+\Delta \tau}} &= {R_\mu}^{t^*}
+ + m_x m_y[\partial_x U'' + \partial_y V'']^{\tau+\Delta \tau}
++ m_y \partial_\eta {\Omega''^{\tau+\Delta \tau}} &= {R_\mu}^{t^*}
\label{mu-small-step}
\\
\delta_\tau \Theta''
-+ m^2[ \partial_x (U''\theta^{t^*})
++ m_x m_y[ \partial_x (U''\theta^{t^*})
+ \partial_y (V''\theta^{t^*})]^{\tau+\Delta \tau}
-+ m \partial_\eta (\Omega''^{\tau+\Delta \tau} \theta^{t^*})
++ m_y \partial_\eta (\Omega''^{\tau+\Delta \tau} \theta^{t^*})
&= {R_\Theta}^{t^*}
\label{theta-small-step}
\\
\delta_\tau W''
-- m^{-1} g\overline{\left[(\alpha/\alpha_d)^{t^*} \biggl[
+- m_y^{-1} g\overline{\left\{(\alpha/\alpha_d)^{t^*} \left[
\partial_\eta (C \partial_\eta \phi'')
- + \partial_\eta\biggl({c_s^2\over\alpha^{t^*}}{\Theta''\over\Theta^{t^*}}\biggr)\biggr]
-- \mu_d''\right]}^\tau
+ +
+\partial_\eta\biggl({c_s^2\over\alpha^{t^*}}{\Theta''\over\Theta^{t^*}}\biggr)
+\right]
+- \mu_d''\right\}}^\tau
&= {R_W}^{t^*}
\label{w-small-step}
\\
\delta_\tau \phi'' + {1\over\mu_d^{t^*}}
-[ m \Omega''^{\tau+\Delta \tau}\phi_\eta^{t^*} - \overline{g W''}^\tau ]
+[ m_y \Omega''^{\tau+\Delta \tau}\phi_\eta^{t^*} - m_y \overline{g W''}^\tau ]
&= {R_\phi}^{t^*}.
\label{geo-small-step}
\end{align}
@@ -180,7 +185,7 @@
%
\begin{align}
R_U^{t^*} = &
-- m[\partial_x(Uu) + \partial_y(Vu)] - \partial_\eta (\Omega u)
+- m_x[\partial_x(Uu) + \partial_y(Vu)] - \hphantom{(m_y/m_x)} \partial_\eta (\Omega u)
- ({\mu}_d \alpha \partial_x p'
- {\mu}_d \alpha' \partial_x \bar{p}) ~~~~~~~~ </font>
<font color="gray">otag
\\
@@ -191,7 +196,7 @@
\\
%
R_V^{t^*} = &
-- m[\partial_x (Uv) + \partial_y (Vv)] - \partial_\eta (\Omega v)
+- m_y[\partial_x (Uv) + \partial_y (Vv)] - (m_y/m_x) \partial_\eta (\Omega v)
- ({\mu}_d \alpha \partial_y p'
- {\mu}_d \alpha' \partial_y \bar{p}) ~~~~~~~~ </font>
<font color="gray">otag
\\
@@ -202,29 +207,29 @@
\\
%
R_{\mu_d}^{t^*} = &
-- m^2[\partial_x U + \partial_y V] - m \partial_\eta \Omega
+- m_x m_y[\partial_x U + \partial_y V] - m_y \partial_\eta \Omega
\label{mu-rhs}
\\
R_\Theta^{t^*} = &
-- m^2[\partial_x (U\theta) + \partial_y (V\theta)] - m \partial_\eta
+- m_x m_y [\partial_x (U\theta) + \partial_y (V\theta)] - m_y \partial_\eta
(\Omega \theta) + F_\Theta
\label{theta-rhs}
\\
%
R_W^{t^*} = &
-- m[\partial_x (Uw) + \partial_y (Vw)] - \partial_\eta
+- (m_x m_y/m_y) [\partial_x (Uw) + \partial_y (Vw)] - \partial_\eta
(\Omega w) ~~~~~~~~ ~~~~~~~~ ~~~~~~~~ ~~~~~~~~ ~~~ </font>
<font color="gray">otag
\\
-& + m^{-1} g (\alpha/\alpha_d) [\partial_\eta p'
+& + m_y^{-1} g (\alpha/\alpha_d) [\partial_\eta p'
- {\bar{\mu}}_d (q_v + q_c +q_r)]
- - m^{-1} {\mu}_d'g + F_W
+ - m_y^{-1} {\mu}_d'g + F_W
\label{w-rhs}
\\
%
R_\phi^{t^*} = &
- \mu_d^{-1}
-[m^2 (U\phi_x + V\phi_y) + m
-\Omega\phi_\eta -gW ],
+[m_x m_y (U\phi_x + V\phi_y) + m_y
+\Omega\phi_\eta - m_y gW ],
\label{geo-rhs}
%
\end{align}
@@ -270,11 +275,12 @@
\Omega''$ term such that
%
\begin{equation}
- \delta_\tau\mu_d = m^2 \int_1^0
+ \delta_\tau\mu_d = m_x m_y \int_1^0
[\partial_x U'' + \partial_y V'']^{\tau+\Delta \tau}
d\eta.
\label{omega}
\end{equation}
+%
After computing ${\mu_d''}^{\tau+\Delta\tau}$ from \eqref{omega},
${\Omega''}^{\tau+\Delta\tau}$ is recovered by using
\eqref{mu-small-step}
@@ -294,21 +300,10 @@
${\alpha_d''}^{\tau+\Delta\tau}$ are recovered from \eqref{p-linear}
and \eqref{small-hydro}.
-\subsection{Full Time-Split Integration Sequence}
-\label{full_time_split_integration}
-
-The time-split RK3 integration technique is summarized below.
-It consists of two primary loops--- an outer loop for the
-large-time-step Runge-Kutta integration, and an inner loop for the
-acoustic mode integration.
-%\bigskip
-
</font>
<font color="red">oindent
\begin{figure}[h]
-%\begin{figure}
\setlength{\fboxrule}{.75pt}
\framebox[\columnwidth]{
-%\parbox{\columnwidth}{
\parbox{6.5truein}{
\vskip 5truept
</font>
<font color="gray">oindent
@@ -344,15 +339,26 @@
(3) Advance horizontal momentum, \eqref{u-small-step}
and \eqref{v-small-step} \hfill\break
%
+\hphantom{BeginBeginBegini(3) }
+Global: Apply polar filter to $U''^{\tau+\Delta \tau}$
+and $V''^{\tau+\Delta \tau}$.
+\hfill\break
+%
\hphantom{BeginBeginBegini}
(4) Advance $\mu_d$ \eqref{mu-small-step} and compute
${\Omega''}^{\tau+\Delta \tau}$ then
advance $\Theta$
\eqref{theta-small-step} \hfill\break
+\hphantom{BeginBeginBegini(4) }
+Global: Apply polar filter to $\mu_d^{\tau+\Delta \tau}$
+and $\Theta''^{\tau+\Delta \tau}$. \hfill\break
%
\hphantom{BeginBeginBegini}
(5) Advance $W$ and $\phi$ \eqref{w-small-step}
and \eqref{geo-small-step} \hfill\break
+\hphantom{BeginBeginBegini(5) }
+Global: Apply polar filter to $W''^{\tau+\Delta \tau}$
+and $\phi''^{\tau+\Delta \tau}$. \hfill\break
%
\hphantom{BeginBeginBegini}
(6) Diagnose $p''$ and $\alpha''$ using
@@ -369,7 +375,9 @@
substep \eqref{rk3a}, \eqref{rk3b} or \eqref{rk3c} \hfill \break
\hphantom{BeginBeginBegin}
(using mass fluxes $U$, $V$ and $\Omega$ time-averaged over the acoustic steps).
-\smallskip \hfill \break
+\hfill \break
+\hphantom{BeginBeginBegin}
+Global: Apply polar filter to scalars. \smallskip\hfill\break
%
\hphantom{BeginBegin}
(8) Compute $p'$ and $\alpha'$ using
@@ -382,15 +390,10 @@
%
\smallskip
\hphantom{Begin}
-%(9) Scalar transport option 2: Advance scalars
-%\eqref{scalar-pert} with a single forward step
-%\hfill \break
-%\hphantom{BeginBegin}
-%(using mass fluxes $U$, $V$ and $\Omega$ time-averaged over the final RK3 acoustic steps).
-%\smallskip \hfill \break
-%\hphantom{Begin}
(9) Compute non-RK3 physics (currently microphysics), advance variables.
-\medskip \hfill \break
+\hfill \break
+\hphantom{Begin(4) }
+Global: Apply polar filter to updated variables. \medskip\hfill\break
%
{\bf End Time Step}
\vskip 5truept
@@ -403,6 +406,14 @@
\label{time_integration_figure}
\end{figure}
+\subsection{Full Time-Split Integration Sequence}
+\label{full_time_split_integration}
+
+The time-split RK3 integration technique is summarized in Figure \ref{time_integration_figure}.
+It consists of two primary loops--- an outer loop for the
+large-time-step Runge-Kutta integration, and an inner loop for the
+acoustic mode integration.
+
In the RK3 scheme, physics can be integrated within the RK3 time
integration (using a time forward step, i.e., step (1) in Fig.
\ref{time_integration_figure}, or the RK3 time integration if higher
@@ -470,20 +481,21 @@
is retained, including the acoustic
step loop. Steps (5) and (6) in the acoustic-step loop, where
$W$ and $\phi$ are advanced and $p''$ and $\alpha''$ are diagnosed, are
-replaced by (1), the diagnosis of the hydrostatic pressure using the definition of
+replaced by the following three steps.
+(1) Diagnose the hydrostatic pressure using the definition of
the vertical coordinate
%
\begin{equation}
-\delta_\eta p_h = {\alpha_d \over \alpha} \mu_d = \bigl(1 + \sum q_m \bigr)\mu_d,
+\delta_\eta p_h = {\alpha_d \over \alpha} \mu_d = \bigl(1 + \sum q_m \bigr)\mu_d.
</font>
<font color="black">otag
\end{equation}
%
</font>
<font color="gray">oindent
-followed by (2), the diagnosis of $\alpha_d$ using
+(2) Diagnose $\alpha_d$ using
the equation of state \eqref{moist-state-equation}
-and the prognosed $\theta$,
-and (3), the diagnosis
-of the geopotential using the hydrostatic equation
+and the prognosed $\theta$.
+(3) Diagnose
+the geopotential using the hydrostatic equation
%
\begin{equation}
\delta_\eta \phi' = - \bar \mu_d \alpha_d' - \mu_d' \alpha_d.
@@ -494,16 +506,16 @@
The vertical velocity
$w$ can be diagnosed from the geopotential equation, but it is not needed
in the solution procedure. The acoustic step loop advances gravity waves,
-including the external mode, when the hydrostatic option is used; there are
-no horizontally propagating acoustic modes in this hydrostatic system.
+including the external mode, and the Lamb wave
+when the hydrostatic option is used.
\section{Spatial Discretization}
The spatial discretization in the ARW solver uses a C grid staggering for
the variables as shown in Fig. \ref{figure:2}. That is, normal
velocities are staggered one-half grid length from the thermodynamic
-variables. The variable indices, $(i,j)$ for the horizontal plane and
-$(i,k)$ for the vertical plane, indicate variable locations where
+variables. The variable indices, $(i,j,k)$
+indicate variable locations with
$(x,y,\eta) = (i\Delta x, j\Delta y, k \Delta \eta)$. We will denote
the points where $\theta$ is located as being {\it mass} points, and
likewise we will denote locations where $u$, $v$, and $w$ are defined as
@@ -525,7 +537,7 @@
we can define the spatial discretization for the ARW solver.
%
-% Figure 3.1
+% Figure 3.something
%
\begin{figure}
\includegraphics *[width=6.5in,bb= 0 0 629.8 325.4]{figures/grids.pdf}
@@ -543,9 +555,9 @@
represented discretely as
%
\begin{equation}
-U = {\mu_d u \over m} \to
-{\overline{\mu_d}^x u \over {\overline m}^x} ,
-~~~ V = {\mu_d v \over m} \to {\overline{\mu_d}^y v \over {\overline m}^y} ,
+U = {\mu_d u \over m_y} \to
+{\overline{\mu_d}^x u \over {\overline m_y}^x} ,
+~~~ V = {\mu_d v \over m_x} \to {\overline{\mu_d}^y v \over {\overline m_x}^y} ,
</font>
<font color="gray">otag
\end{equation}
%
@@ -588,28 +600,30 @@
\label{v-discrete}
\\
\delta_\tau \mu_d''
- + m^2[\delta_x U'' + \delta_y V'']^{\tau+\Delta \tau}
-+ m {\delta_\eta \Omega''^{\tau+\Delta \tau}} & =
+ + m_x m_y[\delta_x U'' + \delta_y V'']^{\tau+\Delta \tau}
++ m_y {\delta_\eta \Omega''^{\tau+\Delta \tau}} & =
R_\mu^{t^*}
\label{mu-discrete}
\\
\delta_\tau \Theta''
-+ m^2[ \delta_x (U''\overline{\theta^{t^*}}^x)
++ m_x m_y[ \delta_x (U''\overline{\theta^{t^*}}^x)
+ \delta_y (V''\overline{\theta^{t^*}}^y)]^{\tau+\Delta \tau}
-+ m \delta_\eta (\Omega''^{\tau+\Delta \tau} \overline{\theta^{t^*}}^\eta)
++ m_y \delta_\eta (\Omega''^{\tau+\Delta \tau} \overline{\theta^{t^*}}^\eta)
&= {R_\Theta}^{t^*}
\label{theta-discrete}
\\
\delta_\tau W''
-- m^{-1} g\overline{\biggl[\overline{(\alpha/\alpha_d)^{t^*}}^\eta
- \delta_\eta (C \delta_\eta \phi'')
- + \delta_\eta\biggl({c_s^2\over\alpha^{t^*}}{\Theta''\over\Theta^{t^*}}\biggr)
-- \mu_d''\biggr]}^\tau
+- m_y^{-1} g\overline{\left\{\overline{(\alpha/\alpha_d)^{t^*}}^\eta
+ \left[ \delta_\eta (C \delta_\eta \phi'')
+ +
+\delta_\eta\biggl({c_s^2\over\alpha^{t^*}}{\Theta''\over\Theta^{t^*}}\biggr)
+\right]
+- \mu_d''\right\}}^\tau
&= {R_W}^{t^*}
\label{w-discrete}
\\
\delta_\tau \phi'' + {1\over\mu_d^{t^*}}
-[m \Omega''^{\tau+\Delta \tau} \delta_\eta \phi^{t^*} - \overline{g W''}^\tau ]
+[m_y \Omega''^{\tau+\Delta \tau} \delta_\eta \phi^{t^*} - m_y \overline{g W''}^\tau ]
&= {R_\phi}^{t^*},
\label{geo-discrete}
\end{align}
@@ -636,7 +650,7 @@
\end{equation}
%
</font>
<font color="gray">oindent
-The operator vertically interpolates variables on mass levels $k$ to the
+This operator vertically interpolates variables on mass levels $k$ to the
$w$ levels $(k+{1\over 2})$. It should be noted that the vertical grid
is defined such that vertical interpolation from $w$ levels to mass
levels reduces to
@@ -655,7 +669,7 @@
- \delta_\eta \overline{\overline{p'}^x}^\eta \delta_x \overline{\phi}^\eta
+ \overline{{\mu}_d'}^x \delta_x \overline{\phi}^\eta ) </font>
<font color="gray">otag \\
& ~~~~~~~~~~~~~~~~~~~~~~~~~ + F_{U_{cor}} + \hbox{advection} +
-\hbox{mixing} + \hbox{physics}
+\hbox{mixing} + \hbox{physics},
\label{u-pg-discrete} \\
%
R_V^{t^*} = &- ( \overline{{\mu}_d}^y \overline{\alpha}^y \delta_y p'
@@ -665,12 +679,12 @@
- \delta_\eta \overline{\overline{p'}^y}^\eta \delta_y \overline{\phi}^\eta
+ \overline{{\mu}_d'}^y \delta_y \overline{\phi}^\eta ) </font>
<font color="red">otag \\
& ~~~~~~~~~~~~~~~~~~~~~~~~~ + F_{V_{cor}} + \hbox{advection} +
-\hbox{mixing} + \hbox{physics}
+\hbox{mixing} + \hbox{physics},
\label{v-pg-discrete} \\
%
-R_W^{t^*} = & ~ m^{-1} g \overline{(\alpha/\alpha_d)}^\eta [\delta_\eta p'
+R_W^{t^*} = & ~ m_y^{-1} g \overline{(\alpha/\alpha_d)}^\eta [\delta_\eta p'
+ {\bar{\mu}}_d \overline{q_m}^\eta]
- - m^{-1} {\mu}_d'g </font>
<font color="blue">otag \\
+ - m_y^{-1} {\mu}_d'g </font>
<font color="gray">otag \\
& ~~~~~~~~~~~~~~~~~~~~~~~~~ + F_{W_{cor}} + \hbox{advection} +
\hbox{mixing} + \hbox{buoyancy} + \hbox{physics}.
\label{w-pg-discrete}
@@ -681,7 +695,8 @@
The terms $F_{U_{cor}}$, $F_{V_{cor}}$, and $F_{W_{cor}}$ in
\eqref{u-pg-discrete} -- \eqref{w-pg-discrete}
represent Coriolis and
-curvature effects in the equations.
+curvature effects in the equations using the isotropic map projections
+(Lambert conformal, polar stereographic, and Mercator) where $m_x=m_y=m$.
These terms in continuous form are given in
\eqref{u-mom-rhs} --
\eqref{w-mom-rhs}. Their spatial discretization is
@@ -697,8 +712,8 @@
{\overline V}^{xy}
- {\overline e}^x
{\overline W}^{x\eta}\, {\overline{\cos \alpha_r}}^x
-- {u {\overline W}^{x\eta} \over r_e}
-</font>
<font color="gray">otag
+- {u {\overline W}^{x\eta} \over r_e},
+\label{coriolis-u-regional}
\\
%
F_{V_{cor}} & = -
@@ -710,8 +725,8 @@
{\overline U}^{xy}
+ {\overline e}^y
{\overline W}^{y\eta}\, {\overline{\sin \alpha_r}}^y
-- {v {\overline W}^{y\eta} \over r_e}
-</font>
<font color="gray">otag
+- {v {\overline W}^{y\eta} \over r_e},
+\label{coriolis-v-regional}
\\
%
F_{W_{cor}} & = + e ({\overline U}^{x\eta}
@@ -720,13 +735,49 @@
{\overline u}^{x\eta} {\overline U}^{x\eta}
+{\overline v}^{y\eta} {\overline V}^{y\eta}
\over r_e}\biggr).
-</font>
<font color="blue">otag
+\label{coriolis-w-regional}
\end{align}
%
</font>
<font color="gray">oindent
Here the operators $\overline{()}^{xy} = \overline{\overline{()}^{x}}^y
$, and likewise for $\overline{()}^{x\eta}$ and $\overline{()}^{y\eta}$.
+For the non-isotropic latitude longitude projection, the Coriolis and
+curvature terms given in
+\eqref{u-mom-rhs-global} and
+\eqref{v-mom-rhs-global} are discretized as
+%
+\begin{align}
+F_{U_{cor}} & = {m_x \over m_y} \biggl[\hphantom{-}{\overline f}^x
+{\overline V}^{xy}
++ {u\overline{V}^{xy} \over r_e} \tan\psi \biggr]
+- {\overline e}^x
+{\overline W}^{x\eta}\, {\overline{\cos \alpha_r}}^x
+- {u {\overline W}^{x\eta} \over r_e},
+\label{coriolis-u-global}
+\\
+%
+F_{V_{cor}} & = {m_y \over m_x} \biggl[ - \overline{f}^y
+\overline{U}^{xy}
+- {\overline{u}^{xy} \overline{U}^{xy} \over r_e} \tan \psi
++ \overline{e}^y \overline{W}^{y\eta}\overline{\sin \alpha_r}^y
++ {v\overline{W}^{y\eta} \over r_e}
+\biggr],
+\label{coriolis-v-global}
+\\
+%
+F_{W_{cor}} & = + e ({\overline U}^{x\eta}
+\cos \alpha_r - (m_x/m_y) {\overline V}^{y\eta}
+\sin \alpha_r) + \biggl({
+ {\overline u}^{x\eta} {\overline U}^{x\eta}
++(m_x/m_y) {\overline v}^{y\eta} {\overline V}^{y\eta}
+\over r_e}\biggr).
+\label{coriolis-w-global}
+%
+\end{align}
+
+
+
\subsection{Advection}
The advection terms in the ARW solver are in the form of a flux divergence and
@@ -735,29 +786,29 @@
%
\begin{align}
R_{U_{adv}}^{t^*} = &
-- m[\partial_x(Uu) + \partial_y(Vu)] + \partial_\eta (\Omega u)
+- m_x[\partial_x(Uu) + \partial_y(Vu)] + \hphantom{(m_x/m_y)} \partial_\eta (\Omega u)
\\
%
R_{V_{adv}}^{t^*} = &
-- m[\partial_x (Uv) + \partial_y (Vv)] + \partial_\eta (\Omega v)
+- m_y[\partial_x (Uv) + \partial_y (Vv)] + (m_x/m_y) \partial_\eta (\Omega v)
\\
%
R_{\mu_{adv}}^{t^*} = &
-- m^2[U_x + V_y] + m \Omega_\eta
+- m_x m_y[U_x + V_y] + m_y \Omega_\eta
\\
R_{\Theta_{adv}}^{t^*} = &
-- m^2[\partial_x (U\theta) + \partial_y (V\theta)] - m \partial_\eta
+- m_x m_y [\partial_x (U\theta) + \partial_y (V\theta)] - m_y \partial_\eta
(\Omega \theta)
\\
%
R_{W_{adv}}^{t^*} = &
-- m[\partial_x (Uw) + \partial_y (Vw)] + \partial_\eta
+- (m_x m_y/m_y) [\partial_x (Uw) + \partial_y (Vw)] + \partial_\eta
(\Omega w)
\\
%
R_{\phi_{adv}}^{t^*} = &
- \mu_d^{-1}
-[m^2 (U\phi_x + V\phi_y) + m
+[m_x m_y (U\phi_x + V\phi_y) + m_y
\Omega\phi_\eta].
\end{align}
%
@@ -767,7 +818,7 @@
%
\begin{equation}
R_{\mu_{adv}}^{t^*} =
-- m^2[\delta_x U + \delta_y V]^{t^*} + m \delta_\eta \Omega^{t^*}.
+- m_x m_y [\delta_x U + \delta_y V]^{t^*} + m_y \delta_\eta \Omega^{t^*}.
\end{equation}
%
@@ -777,12 +828,15 @@
\ref{time_integration_figure}. The spatial discretization used in this
approach is outlined in the next section. For many applications it
is desirable to use positive definite or monotonic advection schemes
-for scalar transport. In the next major release of the ARW we will be
-including a forward-in-time scheme for scalar transport that has
-positive definite and monotonic options. We describe that scheme in
-the section following the description of the RK3 advection.
+for scalar transport.
+%In the next major release of the ARW we will be
+%including a forward-in-time scheme for scalar transport that has
+%positive definite and monotonic options. We describe that scheme in
+%the section following the description of the RK3 advection.
+
\subsubsection{RK3 Advection}
+\label{rk3-advection}
$\hbox{2}^{nd}$ through $\hbox{6}^{th}$ order accurate spatial
discretizations
@@ -796,10 +850,11 @@
%
\begin{equation}
R_{q_{adv}}^{t^*} =
-- m^2[\delta_x (U \overline{q}^{x_{adv}})
+- m_x m_y [\delta_x (U \overline{q}^{x_{adv}})
+ \delta_y (V\overline{q}^{y_{adv}})]
-- m \delta_\eta
+- m_y \delta_\eta
(\Omega \overline{q}^{\eta_{adv}}).
+\label{flux-divergence}
\end{equation}
%
As in the pressure gradient discretization, the discrete operator
@@ -809,6 +864,7 @@
\delta_x (U \overline{q}^{x_{adv}}) = \Delta x^{-1}
\bigl[ (U \overline{q}^{x_{adv}})_{i+1/2} -
(U \overline{q}^{x_{adv}})_{i-1/2} \bigr].
+\label{discrete-divergence}
\end{equation}
%
</font>
<font color="gray">oindent
@@ -866,128 +922,365 @@
The even-order advection operators are spatially centered and thus
contain no implicit diffusion outside of the diffusion inherent in
the RK3 time integration. The odd-order schemes are upwind-biased, and
-the spatial discretization is inherently diffusive. As is evident in
+the spatial discretization is inherently diffusive. The behavior of
+the upwind schemes is easily understood by expanding
+\eqref{discrete-divergence} using the $5^{th}$ order operator, assuming a
+constant mass flux $U$ and
+multiplying by the timestep $\Delta t$:
+%
+\begin{align}
+\Delta t \delta_x (U \overline{q}^{x_{adv}}) =
+& \, {\Delta t} {\delta (Uq)|}^{6th}
+-\left| {U \Delta t \over \Delta x }\right| {1 \over 60}
+\left(-q_{i-3}+6q_{i-2}-15q_{i-1}+20q_{i}-15q_{i+1}+6q_{i+2}-q_{i+3}
+\right)
+</font>
<font color="blue">otag
+\\
+= & \, {\Delta t} {\delta (Uq)|}^{6th} - {Cr \over 60} \Delta x^6 {\partial^6
+q \over \partial x^6} + \, H.O.T.
+</font>
<font color="blue">otag
+\end{align}
+%
+</font>
<font color="blue">oindent
+Similarly, we can expand \eqref{discrete-divergence} using the $3^{rd}$
+order operator:
+\begin{align}
+\Delta t \delta_x (U \overline{q}^{x_{adv}})
+= & \, {\Delta t} {\delta (Uq)|}^{4th} + {Cr \over 12} \Delta x^4 {\partial^4
+q \over \partial x^4} + \, H.O.T.
+</font>
<font color="blue">otag
+\end{align}
+%
+</font>
<font color="red">oindent
+As is evident in
their formulation, the odd-order schemes are comprised of the next
higher (even) order centered scheme plus an upwind term that, for a
constant transport mass flux, is a diffusion term of that next higher
-(even) order with a hyper-viscosity proportional to the Courant number ($Cr$).
+(even) order with a hyper-viscosity proportional to the Courant number
+($Cr$).
+The diffusion term is the leading order error term in the flux
+divergence discretization.
Further details concerning RK3 advection can be found in \citet{wicker02}
-\subsubsection{Forward-In-Time Scalar Advection}
+\subsubsection{Positive-Definite Limiter for RK3 Advection}
+\label{positive-definite-transport}
-A forward-in-time scalar advection scheme having monotonic and positive
-definite options will be available in the next major release of the ARW.
-This option entails bypassing step (7) in the time-split
-integration sequence in Fig. \ref{time_integration_figure}, and adding
-a single advection evaluation before step (9), after the end of the RK3
-loop in the integration sequence. The new advection algorithm is
-patterned after \citet{easter93} and can be written as follows
-%
-\begin{align}
- (\mu_d q)^* &= (\mu_d q)^t + m^2 F_x (q^t)
- \label{split_adv_start} \\
-\mu_d^{*} &= \mu_d^t + m^2 F_x (I) \label{cont-x} \\
-q^{*} &= (\mu_d q)^{*}/ (\mu_d)^{*} \\
-(\mu_d q)^{**} &= (\mu_d q)^* + m^2 F_y (q^*) \\
-\mu_d^{**} &= \mu_d^* + m^2 F_y (I) \label{cont-y} \\
-q^{**} &= (\mu_d q)^{**}/ (\mu_d)^{**} \\
-(\mu_d q)^{t+\Delta t} &= (\mu_d q)^{**} + m F_\eta (q^{**}) \\
-\mu_d^{t+\Delta t} &= \mu_d^{**} + m F_\eta (I) \label{cont-z} \\
-q^{t+\Delta t} & = (\mu_d q)^{t+\Delta t}/\mu_d^{t+\Delta t}. \label{split_adv_end}
-\end{align}
+Mixing ratios of moisture, chemical species or other tracer species
+should remain positive-definite, that is, negative masses should not be
+permitted. The Runge-Kutta transport integration defined by the
+timestepping algorithm \eqref{rk3a} -- \eqref{rk3c}, combined with the
+flux divergence operator \eqref{flux-divergence}, is conservative but it
+does not guarantee positive definiteness; any negative values will be
+offset by positive mass such that mass is conserved. In many physics
+options, negative mixing ratios will be set to zero, and this will
+result in an increase in mass of that species. A positive-definite flux
+renormalization, applied on the final Runge-Kutta transport step
+\eqref{rk3a}, can be used to remove this unphysical effect from the RK3
+scalar transport scheme. This scheme is described in
+\citet{skamarock-weisman-08} and \citet{skamarock2005}, and we briefly
+outline the ARW implementation here.
-</font>
<font color="red">oindent
-In \eqref{split_adv_start} - \eqref{split_adv_end} the operator
+The final RK3 time integration step for transport of a scalar can be
+expressed as
%
\begin{equation}
-F_x (q) = - \Delta t \Delta x^{-1}( (
-{\overline{U}^t q})_{i+1/2} -
-({\overline{U}^t q})_{i-1/2}),
-\label{flux_divergence}
+(\mu\phi)^{t+\Delta t} = (\mu\phi)^{t} - {\Delta t}\bigl\{
+m_x m_y [\delta_x (U \overline{q}^{x_{adv}})
++ \delta_y (V\overline{q}^{y_{adv}})]
++ m_y \delta_\eta
+(\Omega \overline{q}^{\eta_{adv}})\bigr\}
++ \Delta t \mu S_\phi^t
+\label{rk3c-scalar}
\end{equation}
%
</font>
<font color="red">oindent
-with similar definitions
-for $F_y$ and $F_\eta$, and $I$ is a vector with all values equal to
-1.
-The discrete mass continuity
-equation (\eqref{cont-x}+\eqref{cont-y}+\eqref{cont-z}) can be
-written as
+where the flux divergence is evaluated using the (**) time level
+predicted in RK3 step \eqref{rk3b}. The positive-definite flux
+renormalization replaces \eqref{rk3b} with the following two steps. First,
+the scalar mixing ratio is updated using the tendency derived from the
+model physics and source/sink terms.
%
\begin{equation}
-\delta_\tau \mu_d =
-- m^2[\delta_x \overline{U}^t + \delta_y \overline{V}^t]
-- m {\delta_\eta \overline{\Omega}^t}.
-\label{continuity}
+(\mu\phi)^{***} = (\mu\phi)^{t} + {\Delta t} \mu S_\phi^t
+\label{scalar-source-update}
\end{equation}
%
-</font>
<font color="red">oindent
-It can easily be seen that the scheme
-\eqref{split_adv_start} - \eqref{split_adv_end}
-collapses to \eqref{continuity} for $ q = I$,
-and hence it is consistent. The mass fluxes
-$\overline{U}^t$, $\overline{V}^t$, and $\overline{\Omega}^t$,
-represent time-averaged values where the averaging is performed over the
-final RK3 small-time-step cycle. Hence, the mass conservation equation
-\eqref{continuity} produces the mass conservation equation integrated
-within the RK3 scheme. Equation \eqref{continuity} is re-integrated
-within the time-split transport scheme only because a consistent
-column mass $\mu_d$ is
-needed on the sub-steps to retrieve $ q$ from the prognostic
-variable $\mu_d q$.
+where we denote this new predictor as $(\mu\phi)^{***}$. Second, the
+full update is computed using a
+flux divergence composed of a first-order upwind flux plus a
+higher order correction:
-Any forward-in-time flux operator that is stable
-for $Cr \le 1$ can be used to evaluate
-the operators
-$F_x (q)$, $F_y (q)$,
-or $F_\eta (q)$. A scheme based on the Piecewise Parabolic Method
-(PPM) has been implemented in the ARW. It
-is the non-monotonized PPM advection described in \citet{carpenter1990},
-and it provides the fluxes
-$({\overline{U}^t q})_{i\pm {1\over 2}}$
-used
-in \eqref{flux_divergence} in the flux divergence operators.
-Defining $ q_i$ as the control-volume average mixing ratio for cell
-$i$, zone edge values $q^{adv}$ can be defined as
+\begin{align}
+(\mu\phi)^{t+\Delta t} = (\mu\phi)^{***} - {\Delta t}\bigl\{ \,\,\,
+m_x m_y ( \delta_x &\left[ (U \overline{q}^{x_{adv}})^1
++ R(U \overline{q}^{x_{adv}})' \right]
+</font>
<font color="blue">otag
+\\
++\delta_y &\left[(V\overline{q}^{y_{adv}})^1
++ R(V\overline{q}^{y_{adv}})'\right] )
+</font>
<font color="red">otag
+\\
++ m_y \delta_\eta &\left[
+(\Omega \overline{q}^{\eta_{adv}})^1 + R(\Omega \overline{q}^{\eta_{adv}})' \right]
+\,\,\,\, \bigr\}.
+\label{pd-update}
+\end{align}
%
+In \eqref{pd-update}, $()^1$ denotes a first-order upwind flux and $R()'$
+denotes a renormalized higher-order correction flux. The higher-order
+correction flux is the difference between the full RK3 flux and the
+first-order upwind flux, that is,
+%
\begin{equation}
-q^{adv}_{i+{1\over 2}} = \bigl( 7( q_{i+1}+ q_i) -
-( q_{i+2}+ q_{i-1})\bigr)/12.
+(U \overline{q}^{x_{adv}}) = (U \overline{q}^{x_{adv}})^1 + (U
+\overline{q}^{x_{adv}})',
+\label{correction-def}
\end{equation}
%
-</font>
<font color="red">oindent
-The flux through the $i + {1 \over 2}$ face can be written as
+with similar definitions for $(V \overline{q}^{y_{adv}})'$
+and $(\Omega \overline{q}^{\eta_{adv}})'$. The correction flux is then
+renormalized as follows. First, the upwind fluxes are used to perform a
+partial update of the scalar mass.
%
\begin{equation}
-(\overline{U}^t q)_{i+{1 \over 2}} = {\overline{U}^t}_{i+{1 \over 2}} \bigl[
-q^{adv}_{i+{1 \over 2}} - Cr (q^{adv}_{i+{1 \over 2}} - q_{i})
--Cr(1 - Cr)( q^{adv}_{i-{1 \over 2}} - 2 q_{i}
- + q^{adv}_{i+{1\over 2}}) \bigr]
+(\tilde{\mu\phi}) = (\mu\phi)^{***}
+- {\Delta t} \bigl\{
+m_x m_y [\delta_x (U \overline{q}^{x_{adv}})^1
++ \delta_y (V\overline{q}^{y_{adv}})^1]
++ m_y \delta_\eta
+(\Omega \overline{q}^{\eta_{adv}})^1\bigr\}
+</font>
<font color="red">otag
\end{equation}
%
-</font>
<font color="red">oindent
-for $Cr > 0$, and
+This update is positive definite (and monotonic) by design because it is
+the first-order upwind scheme. Next, a prediction of the minimum
+possible values of the new-time-level species mass is computed at each
+point by using only the outward directed fluxes (fluxes that remove mass
+from the control volume),
%
\begin{equation}
-(\overline{U}^t q)_{i+{1 \over 2}} = {\overline{U}^t}_{i+{1 \over 2}} \bigl[
-q^{adv}_{i+{1 \over 2}} - Cr ( q^{adv}_{i+{1 \over 2}}- q_{i+1})
- + Cr(1 + Cr)( q^{adv}_{i+{1 \over 2}} - 2 q_{i+1}
- + q^{adv}_{i+{3\over 2}}) \bigr]
+(\mu\phi)_{min}^{t+\Delta t} = (\tilde{\mu\phi})
+-{\Delta t} \bigl\{
+m_x m_y [\delta_x (U_+ \overline{q}^{x_{adv}})'
++ \delta_y (V_+\overline{q}^{y_{adv}})']
++ m_y \delta_\eta
+(\Omega_+ \overline{q}^{\eta_{adv}})'\bigr\}
+\label{min-prediction}
\end{equation}
%
-</font>
<font color="red">oindent
-for $Cr < 0$, where $Cr = u_{i+{1\over 2}}\Delta t/\Delta
-x$. In addition, to maintain second-order accuracy for the time-split
-forward-in-time advection scheme, the sequence
-\eqref{split_adv_start} -- \eqref{split_adv_end} is reversed
-every other time step; that is, the flux divergence in $\eta$ is
-computed first followed by $y$ and then $x$.
-Finally, the scheme has also been augmented for Courant numbers
-greater than one using the 1D ``semi-Lagrangian'' flux calculation
-described in \citet{linrood96}. More information about this
-scheme and a description of the limiters
-can be found in \citet{skamarock2005}.
+where $U_+$, $V_+$ and $\Omega_+$ indicated the use of fluxes
+out of a control volume only, that is, only those that contribute to
+lowering the scalar mass. As is obvious from \eqref{min-prediction},
+the scalar mass $(\mu\phi)_{min}^{t+\Delta t} < 0$ if
+%
+\begin{equation}
+(\tilde{\mu\phi}) <
+{\Delta t} \bigl\{
+m_x m_y [\delta_x (U_+ \overline{q}^{x_{adv}})'
++ \delta_y (V_+\overline{q}^{y_{adv}})']
++ m_y \delta_\eta
+(\Omega_+ \overline{q}^{\eta_{adv}})'\bigr\}
+</font>
<font color="blue">otag
+\end{equation}
+%
+For each volume where negative mass is indicated by
+\eqref{min-prediction}, the fluxes are renormalized
+such that the outgoing fluxes and mass in the volume are equivalent.
+%
+\begin{equation}
+R\left(U_+\overline{q}^{x_{adv}}\right)' =
+\left(U_+\overline{q}^{x_{adv}}\right)'
+{ \tilde{\mu \phi} \over
+{\Delta t} \bigl\{
+m_x m_y [\delta_x (U_+ \overline{q}^{x_{adv}})'
++ \delta_y (V_+\overline{q}^{y_{adv}})']
++ m_y \delta_\eta
+(\Omega_+ \overline{q}^{\eta_{adv}})'\bigr\} }
+\label{renormalization}
+\end{equation}
+%
+with a similar renormalization applied to the
+$(V \overline{q}^{y_{adv}})'$
+and $(\Omega \overline{q}^{\eta_{adv}})'$.
+If no negative mass is indicated in \eqref{min-prediction}, the
+correction flux is equal to the
+%
+\begin{equation}
+R\left(U_+\overline{q}^{x_{adv}}\right)' =
+\left(U_+\overline{q}^{x_{adv}}\right)',
+\,\,\,
+R\left(V_+\overline{q}^{y_{adv}}\right)' =
+\left(V_+\overline{q}^{y_{adv}}\right)',
+\,\,\,
+R\left(\Omega_+\overline{q}^{\eta_{adv}}\right)' =
+\left(\Omega_+\overline{q}^{\eta_{adv}}\right)'.
+\label{r-ok}
+\end{equation}
+%
+The renormalized fluxes \eqref{renormalization} and \eqref{r-ok}, along
+with the first order fluxes, are then used in the update equation
+\eqref{pd-update}. Note that if no renormalization is needed, the
+scheme \eqref{pd-update} reverts to the standard RK3 update because
+the definition of the correction \eqref{correction-def}.
+%
+% figure 3.2 for polar b.c
+%
+\begin{figure}
+ \includegraphics *[width=6.5in,bb= 0 0 419 183]{figures/polar_bc.pdf}
+ \caption{\label{figure_pole} Latitude-longitude grid structure in the
+pole region. In the ARW formulation, $r_e\Delta\psi = \Delta y/m_y $.}
+\end{figure}
+
+\subsection{Pole Conditions for the Global Latitude-Longitude Grid}
+
+The latitude-longitude grid has a singularity at the two poles where the
+latitude $\psi = \pm 90^\circ$, as illustrated in Figure
+\ref{figure_pole}. By design, no variable is defined at the pole point.
+The area of the control volume face at the pole is zero, thus a flux at
+the pole point is not needed in the solution of any of the prognostic
+variables. For example, in the finite volume discretization of the mass
+conservation equation \eqref{mu-discrete} for $\mu$ in a control volume
+closest to the pole, the meridional gradient ($y$ gradient) of the mass
+flux $\delta_y V''$ will use zero for the pole contribution to this
+term.
+
+The stencils for advection operators higher than 2nd order, described in
+Section \ref{rk3-advection}, cross the poles for flux calculations at
+the control volume faces located $r_e\Delta\psi$ and $2r_e\Delta\psi$
+from the pole point (Figure \ref{figure_pole}, the $V$ flux is indicated
+in red). In the current implementation of the ARW, we reduce the order
+of the flux operator at these faces so that their stencils do not extend
+across the pole point. While this formally reduces the accuracy of the
+scheme, we have not been able to identify any significant degradation in
+the ARW solutions.
+
+Coriolis and curvature terms are computed for the vertical momentum
+equation \eqref{coriolis-w-global} and the horizontal momentum equation
+for $U$ \eqref{coriolis-u-global}. For the $W$ and $U$ points that lie
+$r_e \Delta \psi/2$ from the pole, the stencils for these terms require a
+value of $V$ at the pole point. We set the value of $V$ at the pole
+equal to the value of $V$ at $\Delta \psi$ to evaluate these
+operators. This approximation is also used for the meridional advection
+of $V$, combined with a lowering of the flux-operator order to avoid
+differencing across the pole (as with the flux divergence terms for the
+other prognostic variables).
+
+%\subsubsection{Forward-In-Time Scalar Advection}
+%
+%A forward-in-time scalar advection scheme having monotonic and positive
+%definite options will be available in the next major release of the ARW.
+%This option entails bypassing step (7) in the time-split
+%integration sequence in Fig. \ref{time_integration_figure}, and adding
+%a single advection evaluation before step (9), after the end of the RK3
+%loop in the integration sequence. The new advection algorithm is
+%patterned after \citet{easter93} and can be written as follows
+%
+%\begin{align}
+% (\mu_d q)^* &= (\mu_d q)^t + m_x m_y F_x (q^t)
+% \label{split_adv_start} \\
+%\mu_d^{*} &= \mu_d^t + m_x m_y F_x (I) \label{cont-x} \\
+%q^{*} &= (\mu_d q)^{*}/ (\mu_d)^{*} \\
+%(\mu_d q)^{**} &= (\mu_d q)^* + m_x m_y F_y (q^*) \\
+%\mu_d^{**} &= \mu_d^* + m_x m_y F_y (I) \label{cont-y} \\
+%q^{**} &= (\mu_d q)^{**}/ (\mu_d)^{**} \\
+%(\mu_d q)^{t+\Delta t} &= (\mu_d q)^{**} + m_y F_\eta (q^{**}) \\
+%\mu_d^{t+\Delta t} &= \mu_d^{**} + m_y F_\eta (I) \label{cont-z} \\
+%q^{t+\Delta t} & = (\mu_d q)^{t+\Delta t}/\mu_d^{t+\Delta t}. \label{split_adv_end}
+%\end{align}
+%
+%</font>
<font color="blue">oindent
+%In \eqref{split_adv_start} - \eqref{split_adv_end} the operator
+%
+%\begin{equation}
+%F_x (q) = - \Delta t \Delta x^{-1}( (
+%{\overline{U}^t q})_{i+1/2} -
+%({\overline{U}^t q})_{i-1/2}),
+%\label{flux_divergence}
+%\end{equation}
+%
+%</font>
<font color="blue">oindent
+%with similar definitions
+%for $F_y$ and $F_\eta$, and $I$ is a vector with all values equal to
+%1.
+%The discrete mass continuity
+%equation (\eqref{cont-x}+\eqref{cont-y}+\eqref{cont-z}) can be
+%written as
+%
+%\begin{equation}
+%\delta_\tau \mu_d =
+%- m_x m_y [\delta_x \overline{U}^t + \delta_y \overline{V}^t]
+%- m_y {\delta_\eta \overline{\Omega}^t}.
+%\label{continuity}
+%\end{equation}
+%
+%</font>
<font color="blue">oindent
+%It can easily be seen that the scheme
+%\eqref{split_adv_start} - \eqref{split_adv_end}
+%collapses to \eqref{continuity} for $ q = I$,
+%and hence it is consistent. The mass fluxes
+%$\overline{U}^t$, $\overline{V}^t$, and $\overline{\Omega}^t$,
+%represent time-averaged values where the averaging is performed over the
+%final RK3 small-time-step cycle. Hence, the mass conservation equation
+%\eqref{continuity} produces the mass conservation equation integrated
+%within the RK3 scheme. Equation \eqref{continuity} is re-integrated
+%within the time-split transport scheme only because a consistent
+%column mass $\mu_d$ is
+%needed on the sub-steps to retrieve $ q$ from the prognostic
+%variable $\mu_d q$.
+%
+%Any forward-in-time flux operator that is stable
+%for $Cr \le 1$ can be used to evaluate
+%the operators
+%$F_x (q)$, $F_y (q)$,
+%or $F_\eta (q)$. A scheme based on the Piecewise Parabolic Method
+%(PPM) has been implemented in the ARW. It
+%is the non-monotonized PPM advection described in \citet{carpenter1990},
+%and it provides the fluxes
+%$({\overline{U}^t q})_{i\pm {1\over 2}}$
+%used
+%in \eqref{flux_divergence} in the flux divergence operators.
+%Defining $ q_i$ as the control-volume average mixing ratio for cell
+%$i$, zone edge values $q^{adv}$ can be defined as
+%
+%\begin{equation}
+%q^{adv}_{i+{1\over 2}} = \bigl( 7( q_{i+1}+ q_i) -
+%( q_{i+2}+ q_{i-1})\bigr)/12.
+%\end{equation}
+%
+%</font>
<font color="blue">oindent
+%The flux through the $i + {1 \over 2}$ face can be written as
+%%
+%\begin{equation}
+%(\overline{U}^t q)_{i+{1 \over 2}} = {\overline{U}^t}_{i+{1 \over 2}} \bigl[
+%q^{adv}_{i+{1 \over 2}} - Cr (q^{adv}_{i+{1 \over 2}} - q_{i})
+%-Cr(1 - Cr)( q^{adv}_{i-{1 \over 2}} - 2 q_{i}
+% + q^{adv}_{i+{1\over 2}}) \bigr]
+%\end{equation}
+%%
+%</font>
<font color="blue">oindent
+%for $Cr > 0$, and
+%
+%\begin{equation}
+%(\overline{U}^t q)_{i+{1 \over 2}} = {\overline{U}^t}_{i+{1 \over 2}} \bigl[
+%q^{adv}_{i+{1 \over 2}} - Cr ( q^{adv}_{i+{1 \over 2}}- q_{i+1})
+% + Cr(1 + Cr)( q^{adv}_{i+{1 \over 2}} - 2 q_{i+1}
+% + q^{adv}_{i+{3\over 2}}) \bigr]
+%\end{equation}
+%
+%</font>
<font color="gray">oindent
+%for $Cr < 0$, where $Cr = u_{i+{1\over 2}}\Delta t/\Delta
+%x$. In addition, to maintain second-order accuracy for the time-split
+%forward-in-time advection scheme, the sequence
+%\eqref{split_adv_start} -- \eqref{split_adv_end} is reversed
+%every other time step; that is, the flux divergence in $\eta$ is
+%computed first followed by $y$ and then $x$.
+%Finally, the scheme has also been augmented for Courant numbers
+%greater than one using the 1D ``semi-Lagrangian'' flux calculation
+%described in \citet{linrood96}. More information about this
+%scheme and a description of the limiters
+%can be found in \citet{skamarock2005}.
+
\section{Stability Constraints}
There are two time steps that a user must specify when running the ARW:
@@ -1003,6 +1296,7 @@
for applications.
\subsection{RK3 Time Step Constraint}
+\label{rk3_timestep_constraint}
The RK3 time step is limited by the advective Courant number $u \Delta
t/\Delta x$ and the user's choice of advection schemes--- users can
@@ -1048,7 +1342,9 @@
in the simulation. For example in real-data applications, where jet
stream winds may reach as high as $100\,\,\hbox{ms}^{-1}$, the maximum
time step would be approximately 80 s on a $\Delta x = 10 \,\,\hbox{km}$ grid
-using $5^{th}$ order advection.
+using $5^{th}$ order advection. For convection-permitting resolutions
+(typically $\Delta x \le 5$ km), the vertical velocities in convective
+updrafts produce the stability-limiting Courant numbers.
Given additional constraint from the time splitting, and to provide a
safety buffer, we usually choose a time step that is approximately 25\%
less than that given by \eqref{rk3_timestep}. This time step is
@@ -1057,7 +1353,7 @@
for choosing a time step is that the time step, in seconds, should be
approximately 3 times the horizontal grid distance, in kilometers.
For the ARW, the time step (in seconds) should be approximately
-6 times the grid distance (in kilometers).
+6 times the grid distance (in kilometers).
\subsection{Acoustic Time Step Constraint}
\label{acoustic_step_constraint}
@@ -1087,3 +1383,102 @@
\eqref{acoustic_timestep}). Note that it is the ratio of the RK3
time step to the acoustic time step that is the required input in the ARW.
+
+\subsection{Adaptive Time Step}
+\label{adaptive_time_step}
+
+The ARW model is typically integrated with a fixed timestep,
+that is chosen to produce a stable integration. During any time in
+the integration, the maximum stable timestep is likely to be larger than
+the fixed timestep. In ARWV3, an adaptive timestepping capability has
+been introduced that chooses the RK3 timestep based on the temporally-evolving
+wind fields. The adaptively-chosen timestep is usually larger than the
+typical fixed timestep, hence the dynamics integrates faster and
+physics are called less often, and the time-to-completion of the
+simulation can be substantially reduced.
+
+In the adaptive timestep scheme, a target maximum Courant number
+$Cr_{target}$ is chosen, where typically $1.1 \leq Cr_{target} \leq
+1.2$. The maximum Courant number in the domain at a given time ($Cr_{domain}$), computed
+for all the velocity components $(u,v,w)$, is then used to compute a new
+timestep. When the maximum Courant number in the domain is less than
+the target maximum Courant number $(Cr_{domain} < Cr_{target})$,
+then the timestep can be increased
+and the new timestep is computed using
+%
+\begin{equation}
+\Delta t_{current} = min \Big( 1+f_{i} ,
+{Cr_{target} \over Cr_{domain}}\Big) \cdot \Delta t_{previous},
+\label{adapt_inc_timestep}
+\end{equation}
+%
+</font>
<font color="blue">oindent
+where a typical value for the regulated increase is $f_{i} \leq 5\%$. When
+the computed maximum domain-wide Courant number exceeds the targeted
+maximum allowable Courant number $(Cr_{domain} > Cr_{target})$,
+then the time step is
+decreased to insure model stability:
+%
+\begin{equation}
+\Delta t_{current} = max \Big( 1-f_{d} ,
+{{Cr_{target} - {0.5 ({Cr_{domain} - Cr_{target}}) }}
+\over Cr_{domain} } \Big) \cdot \Delta t_{previous},
+\label{adapt_dec_timestep}
+\end{equation}
+%
+</font>
<font color="blue">oindent
+where typically the factor to decrease the time step $f_{d} = 25\%$. Both a lower bound and
+an upper bound on the time step are enforced based on the initial settings
+of the time step suggested in section \ref{rk3_timestep_constraint}:
+%
+\begin{equation}
+\Delta t_{init} = 6 \cdot \Delta x ,
+\label{adapt_init_timestep}
+\end{equation}
+\begin{equation}
+\Delta t_{min} = 0.5 \cdot \Delta t_{init},
+\label{adapt_bound_min_timestep}
+\end{equation}
+\begin{equation}
+\Delta t_{max} = 3.0 \cdot \Delta t_{init},
+\label{adapt_bound_max_timestep}
+\end{equation}
+\begin{equation}
+\Delta t = min ( max ( \Delta t, \Delta t_{min} ) , \Delta t_{max} ),
+\label{adapt_bound_min_timestep}
+\end{equation}
+%
+</font>
<font color="gray">oindent
+to guard the time step from becoming too small or too large.
+The computation of the number of acoustic time steps $n_a$, described in section
+\ref{acoustic_step_constraint}, is handled during the first of each of the
+RK3 model integration steps:
+%
+\begin{equation}
+n_a
+= max \Big( 2 \cdot \lfloor { 300 \cdot {\Delta t \over \Delta x} + 1 } \rfloor , 4 \Big).
+\label{adapt_acoustic_timestep}
+\end{equation}
+%
+
+For a simulation using nests, the fine-grid domain must maintain an integer number of
+time steps within a single model integration step from the parent. At each of
+the starting
+model integration steps for the child grid, those integration steps when the time on the fine-grid
+domain equals the time
+on the parent's domain, the adaptive time step algorithm is conducted for the fine grid.
+Note that the ratio of the nominal grid distance between the parent and the child does not
+necesarily imply that same ratio between the model integration time steps.
+
+
+\subsection{Map Projection Considerations}
+\label{global_considerations}
+
+For ARW configurations using the Lambert conformal, polar stereographic,
+or Mercator projections, the timestep constraints is determined by
+the smallest physical horizontal grid spacing,
+i.e. $\hbox{min}(\Delta x/m_x, \Delta y/m_y)$. For global applications,
+the grid distance used to determine the timestep should be $\Delta
+x/m_x$ evaluated at the computational latitude at which the polar
+filters are activated. Polar filtering is discussed in section
+\ref{polar_filter_section}.
Modified: trunk/wrf/technote/equation.tex
===================================================================
--- trunk/wrf/technote/equation.tex 2008-04-26 00:42:07 UTC (rev 61)
+++ trunk/wrf/technote/equation.tex 2008-04-30 17:47:21 UTC (rev 62)
@@ -150,7 +150,8 @@
flux form but we find no advantage in doing so since $\mu \phi$ is not a
conserved quantity. We could also use a prognostic pressure equation in
place of \eqref{prognostic_end} \citep[see][]{laprise92}, but pressure is
-not a conserved variable and we could not use a pressure equation and
+not a conserved variable and we could not use a pressure equation
+together with
the conservation equation for $\Theta$ \eqref{cartesian_theta} because
they are linearly dependent.
Additionally, prognostic pressure
@@ -195,24 +196,24 @@
With these definitions, the moist Euler equations can be written as
%
\begin{align}
-\partial_t U + (</font>
<font color="blue">abla \cdot {\bf V} u)_\eta
+\partial_t U + (</font>
<font color="red">abla \cdot {\bf V} u)
+ \mu_d \alpha \partial_x p
+ (\alpha/\alpha_d) \partial_\eta p \partial_x \phi & = F_U \\
%
-\partial_t V + (</font>
<font color="blue">abla \cdot {\bf V} v )_\eta
+\partial_t V + (</font>
<font color="red">abla \cdot {\bf V} v )
+ \mu_d \alpha \partial_y p
+ (\alpha/\alpha_d) \partial_\eta p \partial_y \phi & = F_V \\
%
-\partial_t W + (</font>
<font color="blue">abla \cdot {\bf V} w )_\eta
+\partial_t W + (</font>
<font color="red">abla \cdot {\bf V} w )
- g [ (\alpha/\alpha_d) \partial_\eta p - \mu_d ] & = F_W \\
%
-\partial_t \Theta + (</font>
<font color="blue">abla \cdot {\bf V} \theta )_\eta & = F_\Theta \\
+\partial_t \Theta + (</font>
<font color="red">abla \cdot {\bf V} \theta ) & = F_\Theta \\
%
-\partial_t \mu_d + (</font>
<font color="blue">abla \cdot {\bf V} )_\eta & = 0 \\
+\partial_t \mu_d + (</font>
<font color="red">abla \cdot {\bf V} ) & = 0 \\
%
-\partial_t \phi + \mu_d^{-1} [({\bf V} \cdot </font>
<font color="blue">abla \phi )_\eta - gW ] & = 0\\
+\partial_t \phi + \mu_d^{-1} [({\bf V} \cdot </font>
<font color="red">abla \phi ) - gW ] & = 0\\
%
-\partial_t Q_m + ({\bf V} \cdot </font>
<font color="blue">abla q_m )_\eta & = F_{Q_m}
+\partial_t Q_m + (</font>
<font color="black">abla \cdot {\bf V} q_m ) & = F_{Q_m}
%
\end{align}
</font>
<font color="gray">oindent
@@ -246,25 +247,31 @@
\section{Map Projections, Coriolis and Curvature Terms}
\label{spherical_projections}
-The ARW solver currently supports three projections to the sphere--- the
-Lambert conformal, polar stereographic, and Mercator projections. These
-projections are described in \citet{haltiner_and_williams}. These
-projections, and the ARW implementation of the map factors, assume that
-the map factor transformations for $x$ are $y$ are identical at a
-given point; that is, the transformation is isotropic. Anisotropic
-transformations, such as a latitude-longitude grid, can be accommodated
-by defining separate map factors for the $x$ and $y$ transformations.
+The ARW solver currently supports four projections to the sphere--- the
+Lambert conformal, polar stereographic, Mercator, and latitude-longitude
+projections. These projections are described in
+\citet{haltiner_and_williams}. The transformation is isotropic for
+three of these projections -- the Lambert conformal, polar stereographic,
+and Mercator grids. An isotropic transformation requires $(\Delta
+x/\Delta y)|_{earth} = \hbox{constant}$ everywhere on the grid. Only
+isotropic transformations were supported in the previous ARW releases.
+Starting with the ARWV3 release, we now support anisotropic
+projections, in this case the latitude-longitude grid, and with it the
+full latitude-longitude global model. The ARW implements the
+projections using map factors, and the generalization to anisotropic
+transformations introduced in ARW V3 requires that there be map factors
+for both the $x$ and $y$ components of the transformation from
+computational to physical space in order to accomodate the anisotropy.
In the ARW's computational space, $\Delta x$ and $\Delta y$ are
constants. Orthogonal projections to the sphere require that
the physical distances between grid points in the projection vary with
-position in the grid. To transform the governing equations, a
-map scale factor $m$ is defined as the ratio of the distance
-in computational space to the
-corresponding distance on the earth's surface:
+position on the grid. To transform the governing equations,
+map scale factors $m_x$ and $m_y$ are defined as the ratio of the distance
+in computational space to the corresponding distance on the earth's surface:
\begin{equation}
-m = {(\Delta x, \Delta y) \over \hbox{distance on the earth}}.
+(m_x,m_y) = {(\Delta x, \Delta y) \over \hbox{distance on the earth}}.
\end{equation}
%
</font>
<font color="gray">oindent
@@ -272,7 +279,8 @@
by redefining the momentum variables as
%
\begin{equation}
-U = {\mu}_d u/m,~~~V = {\mu}_d v/m, ~~~W = \mu_d w / m, ~~~\Omega = \mu_d \dot\eta /m.
+U = {\mu}_d u/m_y,~~~V = {\mu}_d v/m_x, ~~~W = \mu_d w / m_y,
+~~~\Omega = \mu_d \dot\eta /m_y.
</font>
<font color="gray">otag
\end{equation}
%
@@ -282,39 +290,41 @@
can be written as
%
\begin{align}
-\partial_t U + m[\partial_x(Uu) + \partial_y(Vu)] + \partial_\eta (\Omega u)
+\partial_t U + m_x[\partial_x(Uu) + \partial_y(Vu)] +
+\hphantom{(m_y/m_x)}
+\partial_\eta (\Omega u)
+ \mu_d \alpha \partial_x p
+ (\alpha/\alpha_d) \partial_\eta p \partial_x \phi & = F_U
\label{u-mom-full}
\\
%
-\partial_t V + m[\partial_x (Uv) + \partial_y (Vv)] + \partial_\eta (\Omega v)
+\partial_t V + m_y[\partial_x (Uv) + \partial_y (Vv)] + (m_y/m_x)\partial_\eta (\Omega v)
+ \mu_d \alpha \partial_y p
+ (\alpha/\alpha_d) \partial_\eta p \partial_y \phi & = F_V
\label{v-mom-full}
\\
%
-\partial_t W + m[\partial_x (Uw) + \partial_y (Vw)] + \partial_\eta (\Omega w)
-- m^{-1} g [ (\alpha/\alpha_d) \partial_\eta p - \mu_d ] & = F_W
+\partial_t W + (m_x m_y/m_y) [\partial_x (Uw) + \partial_y (Vw)] + \partial_\eta (\Omega w)
+- m_y^{-1} g [ (\alpha/\alpha_d) \partial_\eta p - \mu_d ] & = F_W
\label{w-mom-full}
\\
%
\partial_t \Theta +
-m^2[\partial_x (U\theta) + \partial_y (V\theta)] + m \partial_\eta (\Omega \theta)
+m_x m_y[\partial_x (U\theta) + \partial_y (V\theta)] + m_y \partial_\eta (\Omega \theta)
& = F_\Theta \\
%
-\partial_t \mu_d + m^2[U_x + V_y] + m \partial_\eta (\Omega)
+\partial_t \mu_d + m_x m_y[U_x + V_y] + m_y \partial_\eta (\Omega)
& = 0
\label{mass-cons-full}
\\
%
-\partial_t \phi + \mu_d^{-1} [m^2 (U\phi_x + V\phi_y) + m
-\Omega\phi_\eta - gW ] & = 0
+\partial_t \phi + \mu_d^{-1} [m_x m_y (U\phi_x + V\phi_y) + m_y
+\Omega\phi_\eta - m_y gW ] & = 0
\label{geo-full}
\\
%
\partial_t Q_m +
-m^2[\partial_x (U q_m) + \partial_y (V q_m)] + m \partial_\eta (\Omega q_m)
+m_x m_y \partial_x (U q_m) + \partial_y (V q_m)] + m _y\partial_\eta (\Omega q_m)
& = F_{Q_m},
%
\end{align}
@@ -339,8 +349,10 @@
The right-hand-side terms of the momentum equations
\eqref{u-mom-full} -- \eqref{w-mom-full} contain the Coriolis and
curvature terms along with mixing terms and
-physical forcings. The Coriolis and curvature terms
-can be written as follows:
+physical forcings. For the isotropic projections
+(Lambert conformal, polar stereographic, and Mercator), where
+$m_x=m_y=m$, the Coriolis and
+curvature terms are cast in the following form:
\begin{align}
F_{U_{cor}} & = + \biggl(f + u {\partial m \over
\partial y} - v {\partial m \over \partial x}\biggr) V
@@ -350,9 +362,10 @@
%
F_{V_{cor}} & = - \biggl(f + u {\partial m \over
\partial y} - v {\partial m \over \partial x}\biggr) U
-+ e W \sin \alpha_r - {vW \over r_e} \\
++ e W \sin \alpha_r - {v W \over r_e} \\
%
-F_{W_{cor}} & = + e (U \cos \alpha_r - V \sin \alpha_r) + \biggl({ uU +vV \over r_e}\biggr),
+F_{W_{cor}} & = + e (U \cos \alpha_r - V \sin \alpha_r)
++ \biggl({ uU +vV \over r_e}\biggr),
\label{w-mom-rhs}
\end{align}
</font>
<font color="gray">oindent
@@ -367,15 +380,37 @@
containing $r_e$ relate to vertical (earth-surface) curvature, and those
with $e$ and $f$ are the Coriolis force.
-In idealized cases, the map scale factor $m = 1$,
+The curvature and Coriolis terms for the momentum equations
+are cast in the following form
+for the anisotropic latitude-longitude grid:
+%
+\begin{align}
+F_{U_{cor}} & = {m_x \over m_y} \biggl[fV + {uV \over r_e}
+\tan\psi \biggr] - {uW \over r_e} -eW\cos \alpha_r
+\label{u-mom-rhs-global}
+\\
+%
+F_{V_{cor}} & = {m_y \over m_x} \biggl[ - fU - {uU \over r_e}
+\tan \psi + {vW \over r_e}
++ eW\sin \alpha_r \biggr]
+\label{v-mom-rhs-global}
+\\
+%
+F_{W_{cor}} & = + e (U \cos \alpha_r - (m_x/m_y) V \sin \alpha_r)
++ \biggl({ uU + (m_x/m_y) vV \over r_e}\biggr),
+\label{w-mom-rhs-global}
+%
+\end{align}
+
+For idealized cases, the map scale factor $m_x = m_y = 1$,
$f$ is often taken to be constant, and $e = 0$.
\section{Perturbation Form of the Governing Equations}
Before constructing the discrete solver, it is advantageous to recast
-the governing equations using perturbation variables so as to reduce
+the governing equations using perturbation variables to reduce
truncation errors in the horizontal pressure gradient calculations in
-the discrete solver, in addition to reducing machine rounding errors in
+the discrete solver and machine rounding errors in
the vertical pressure gradient and buoyancy calculations.
For this purpose, new
variables are defined as perturbations from a hydrostatically-balanced
@@ -383,8 +418,9 @@
variables (denoted by overbars) that are a function of height only and
that satisfy the governing equations for an atmosphere at rest. That is,
the reference state is in hydrostatic balance and is strictly only a
-function of $z$. In this manner, $p=\bar p(z)+p'$, $\phi=\bar \phi(z)
-+\phi'$, $\alpha=\bar \alpha(z) +\alpha'$, and $\mu_d = \bar\mu_d(x,y) +
+function of $\bar z$. In this manner, $p=\bar p(\bar z)+p'$, $\phi=\bar
+\phi(\bar z)
++\phi'$, $\alpha=\bar \alpha(\bar z) +\alpha'$, and $\mu_d = \bar\mu_d(x,y) +
\mu_d'$. Because the $\eta$ coordinate surfaces are generally not
horizontal, the reference profiles $\bar p$, $\bar\phi$, and
$\bar\alpha$ are functions of $(x,y,\eta)$.
@@ -395,7 +431,9 @@
\eqref{u-mom-full} -- \eqref{w-mom-full} are written as
%
\begin{align}
- \partial_t U + m[\partial_x(Uu) + \partial_y(Vu)] + \partial_\eta (\Omega u)
+ \partial_t U + m_x[\partial_x(Uu) + \partial_y(Vu)] +
+\hphantom{(m_y/m_x)}
+\partial_\eta (\Omega u)
+ ({\mu}_d \alpha \partial_x p'
+ {\mu}_d \alpha' \partial_x \bar{p}) ~~~~~~~~ </font>
<font color="gray">otag &
\\
@@ -406,7 +444,7 @@
\label{u-mom-pert}
\\
%
-\partial_t V + m[\partial_x (Uv) + \partial_y (Vv)] + \partial_\eta (\Omega v)
+\partial_t V + m_y[\partial_x (Uv) + \partial_y (Vv)] + (m_y/m_x) \partial_\eta (\Omega v)
+ ({\mu}_d \alpha \partial_y p'
+ {\mu}_d \alpha' \partial_y \bar{p}) ~~~~~~~~ </font>
<font color="gray">otag &
\\
@@ -415,12 +453,12 @@
- {\mu}_d' \partial_y \phi )
&= F_V \\
%
-\partial_t W + m[\partial_x (Uw) + \partial_y (Vw)] + \partial_\eta
+\partial_t W + (m_x m_y/m_y) [\partial_x (Uw) + \partial_y (Vw)] + \partial_\eta
(\Omega w) ~~~~~~~~ ~~~~~~~~ ~~~~~~~~ ~~~~~~~~ ~~~ </font>
<font color="gray">otag &
\\
-- m^{-1} g (\alpha/\alpha_d) [\partial_\eta p'
+- m_y^{-1} g (\alpha/\alpha_d) [\partial_\eta p'
- {\bar{\mu}}_d (q_v + q_c +q_r)]
-+ m^{-1} {\mu}_d'g
++ m_y^{-1} {\mu}_d'g
&= F_W,
\label{w-mom-pert}
\end{align}
@@ -429,14 +467,14 @@
and geopotential equation \eqref{geo-full} become
%
\begin{align}
-\partial_t \mu_d' + m^2[\partial_x U + \partial_y V] + m
+\partial_t \mu_d' + m_x m_y[\partial_x U + \partial_y V] + m_y
\partial_\eta \Omega
& = 0 \\
%
\partial_t \phi'
+ \mu_d^{-1}
-[m^2 (U\phi_x + V\phi_y) + m
-\Omega\phi_\eta -gW ] = 0.
+[m_x m_y (U\phi_x + V\phi_y) + m_y
+\Omega\phi_\eta -m_y gW ] = 0.
%
\end{align}
%
@@ -446,11 +484,11 @@
scalars
\begin{align}
\partial_t \Theta
-+ m^2[\partial_x (U\theta) + \partial_y (V\theta)] + m \partial_\eta (\Omega \theta)
++ m_x m_y[\partial_x (U\theta) + \partial_y (V\theta)] + m_y \partial_\eta (\Omega \theta)
& = F_\Theta \\
%
\partial_t Q_m
-+ m^2[\partial_x (U q_m) + \partial_y (V q_m)] + m \partial_\eta (\Omega q_m)
++ m_x m_y[\partial_x (U q_m) + \partial_y (V q_m)] + m_y \partial_\eta (\Omega q_m)
& = F_{Q_m}.
\label{scalar-pert}
\end{align}
Modified: trunk/wrf/technote/filter.tex
===================================================================
--- trunk/wrf/technote/filter.tex 2008-04-26 00:42:07 UTC (rev 61)
+++ trunk/wrf/technote/filter.tex 2008-04-30 17:47:21 UTC (rev 62)
@@ -4,7 +4,9 @@
A number of formulations for turbulent mixing and filtering are
available in the ARW solver. Some of these filters are used for
numerical reasons. For example, divergence damping is used to filter
-acoustic modes from the solution. Other filters are meant to represent
+acoustic modes from the solution and polar filtering is used to reduce
+the timestep restriction arising from the converging gridlines of the
+latitude-longitude grid. Other filters are meant to represent
sub-grid turbulence processes that cannot be resolved on the chosen
grid. These filters remove energy from the solution and are formulated
in part on turbulence theory and observations, or represent energy sink
@@ -21,15 +23,91 @@
in this chapter, other numerical filters available in the ARW solver are
described.
+
+\section{Latitude-Longitude Global Grid and Polar Filtering}
+\label{polar_filter_section}
+
+Polar filtering is used to reduce the timestep restriction associated
+with the gridlines that converge as they approach the poles of the
+latitude-longitude grid; the converging gridlines reduce the
+longitudinal gridlength, and thus the stable timestep for transport,
+discussed in Section \ref{rk3_timestep_constraint}, and for acoustic
+waves, discussed in Section \ref{acoustic_step_constraint}, will be reduced.
+
+Polar filtering a given variable is accomplished by first applying a 1D
+Fourier transform to the variable on a constant computational-grid
+latitude circle (the forward transform), where the field is periodic.
+The Fourier coefficients with wavenumbers above a prescribed threshold
+are truncated, after which a transformation back to physical space is
+applied (the backward transform), completing the filter step.
+
+A filter application can be written as
+%
+\begin{equation}
+\hat\phi(k)_{filtered} = a(k)\,\hat\phi(k), \,\,\,\,\, \hbox{for all }k,
+</font>
<font color="blue">otag
+\end{equation}
+%
+where $\hat\phi(k)$ and $\hat\phi(k)_{filtered}$ are the Fourier
+coefficients for the generic variable $\phi$ before and after filtering,
+and $a(k)$ are the filter coefficients defined as a function of the
+dimensionless wavenumber $k$. The ARW polar filter coefficients $a(k)$
+for the Fourier amplitudes are
+%
+\begin{equation}
+a(k) = \hbox{min}\left[1.,\hbox{max}\left(0.,
+\left({\cos{\psi} \over \cos{\psi_o}}\right)^2 {1 \over \sin^2({\pi k/n)}}
+\right)\right]
+</font>
<font color="gray">otag
+%\label{polar-filter}
+\end{equation}
+%
+where $\psi$ is the latitude on the computational grid, $\psi_o$ is the
+latitude above which the polar filter is applied (the filter is applied
+for $ |\psi| > \psi_o$), and $n$ is the number of grid points in the latitude
+circle. For even $n$, the grid
+admits wavenumber $k=0$ and dimensionless wavenumbers $\pm k$ for
+$k=1\to n/2 - 1$, and the $2\Delta$ wave $k=n/2$. For odd $n$, the grid
+admits wavenumber $k=0$, and dimensionless wavenumbers $\pm k$ for
+$k=1\to (n-1)/2 $.
+
+Polar filter applications within the split-explicit RK3 time integration
+scheme are given in Figure \ref{time_integration_figure}.
+Within the acoustic steps,
+the filter is applied to
+${U''}^{\tau+\Delta\tau}$ and ${V''}^{\tau+\Delta\tau}$ immediately
+after the horizontal momentum is advanced, to $\mu_d^{\tau+\Delta\tau}$
+and ${\Theta''}^{\tau+\Delta\tau}$ immediately after the column mass and
+the potential temperature are advanced, and to ${W''}^{\tau+\Delta\tau}$
+and ${\phi''}^{\tau+\Delta\tau}$ after the vertically implicit part of
+the acoustic step is complete. The polar filter is also applied to the
+scalars after they are advanced every RK3 sub-step
+\eqref{rk3a}-\eqref{rk3c} and after the microphysics step.
+
+The prognostic variables are coupled with the column mass $\mu$ when
+filtered, and in this way mass is conserved. An exception to this is
+the geopotential $\phi$ which is filtered without being coupled. The
+Fourier filtering is conservative, but it is not monotonic or positive
+definite. As noted in the ARW stability discussion in Section
+\ref{global_considerations}, timestep stability should be based on the longitudinal
+gridlength where the polar filters are activated, $(\Delta
+x/m_x)|_{\psi_o}$. The positive definite transport option, presented in
+Section \ref{positive-definite-transport}, will not necessarily produce
+positive-definite results because of the converging gridlines and should
+not be used.
+
\section{Explicit Spatial Diffusion}
-The ARW solver has two formulations for spatial dissipation--- diffusion
-along coordinate surfaces and diffusion in physical $(x,y,z)$ space. In
-the following sections we present the diffusion operators for these two
+The ARW solver has three formulations for spatial dissipation--- diffusion
+along coordinate surfaces, diffusion in physical $(x,y,z)$ space, and
+a sixth-order diffusion applied on horizontal coordinate surfaces. In
+the following sections we present the diffusion operators for the the
+first two
formulations, followed by the four separate formulations that can be
-used to compute the eddy viscosities. We conclude with a description of
+used to compute the eddy viscosities and a description of
the prognostic turbulent kinetic energy (TKE) equation used in one set
-of these formulations.
+of these formulations. The sixth order spatial filter is described at
+the end of this section.
\subsection{Horizontal and Vertical Diffusion on Coordinate Surfaces}
@@ -42,8 +120,8 @@
%
\begin{equation}
\partial_t (\mu_d a) = ... \,\,
-+ \mu_d \bigl[ m \partial_x(m K_h \partial_x a)
- + m \partial_y(m K_h \partial_y a)\bigr]
++ \mu_d \bigl[ m_x \partial_x(m_x K_h \partial_x a)
+ + m_y \partial_y(m_y K_h \partial_y a)\bigr]
+ g^2(\mu_d\alpha)^{-1} \partial_\eta(K_v \alpha^{-1} \partial_\eta a).
\label{dissipation_1}
\end{equation}
@@ -53,30 +131,30 @@
%
\begin{align}
\partial_t U &= ... \,\, +
-{\overline\mu_d}^x {\overline m}^x
-\bigl[ \delta_x(m K_h \delta_x u)
-+ \delta_y({\overline m}^{xy} {\overline K_h}^{xy} \delta_y u) \bigr]
-+ g^2({\overline \mu_d}^x {\overline\alpha}^x)^{-1}
+{\overline{\mu_d}}^x {\overline{m_x}}^x
+\bigl[ \delta_x(m_x K_h \delta_x u)
++ \delta_y({\overline{m_y}}^{xy} {\overline{K_h}}^{xy} \delta_y u) \bigr]
++ m_y^{-1} g^2({\overline{\mu_d}}^x {\overline\alpha}^x)^{-1}
\delta_\eta(K_v ({\overline \alpha}^{x\eta})^{-1} \delta_\eta u)
</font>
<font color="red">otag \\
%
\partial_t V &= ... \,\, +
-{\overline\mu_d}^y {\overline m}^y
+{\overline{\mu_d}}^y {\overline{m_y}}^y
\bigl[
-\delta_x({\overline m}^{xy} {\overline K_h}^{xy} \delta_x v)
-+\delta_y(m K_h \delta_y v)
+\delta_x({\overline{m_x}}^{xy} {\overline K_h}^{xy} \delta_x v)
++\delta_y(m_y K_h \delta_y v)
\bigr]
-+ g^2({\overline \mu_d}^y {\overline\alpha}^y)^{-1}
++ m_x^{-1} g^2({\overline{\mu_d}}^y {\overline\alpha}^y)^{-1}
\delta_\eta(K_v ({\overline \alpha}^{y\eta})^{-1} \delta_\eta v)
</font>
<font color="red">otag \\
%
\partial_t W &= ... \,\, +
-\mu_d m
+\mu_d m_x
\bigl[
-\delta_x({\overline m}^{x} {\overline K_h}^{x\eta} \delta_x w)
-+\delta_y({\overline m}^{y} {\overline K_h}^{y\eta} \delta_y w)
+\delta_x({\overline m_x}^{x} {\overline K_h}^{x\eta} \delta_x w)
++\delta_y({\overline m_y}^{y} {\overline K_h}^{y\eta} \delta_y w)
\bigr]
-+ g^2({\mu_d} {\overline\alpha}^\eta)^{-1}
++ m_y^{-1} g^2({\mu_d} {\overline\alpha}^\eta)^{-1}
\delta_\eta(K_v \alpha^{-1} \delta_\eta w).
</font>
<font color="gray">otag
\end{align}
@@ -87,10 +165,10 @@
%
\begin{equation}
\partial_t (\mu_d q) = ... \,\, +
-\mu_d m
+\mu_d m_x m_y
\bigl[
-\delta_x({\overline m}^{x} P_r^{-1} {\overline K_h}^{x} \delta_x q)
-+\delta_y({\overline m}^{y} P_r^{-1} {\overline K_h}^{y} \delta_y q)
+\delta_x({\overline m_x}^{x} P_r^{-1} {\overline K_h}^{x} \delta_x q)
++\delta_y({\overline m_y}^{y} P_r^{-1} {\overline K_h}^{y} \delta_y q)
\bigr]
+ g^2({\mu_d} {\alpha})^{-1}
\delta_\eta(K_v \alpha^{-1} \delta_\eta q).
@@ -138,19 +216,19 @@
%
\begin{align}
\partial_t U &= ...\,\,
-- m \bigl[ \partial_x \tau_{11} + \partial_y \tau_{12}
+- m_x \bigl[ \partial_x \tau_{11} + \partial_y \tau_{12}
- \partial_z ( z_x \tau_{11} + z_y \tau_{12}) \bigr]
- \partial_z \tau_{13}
\label{u-disp-phys} \\
%
\partial_t V &= ...\,\,
-- m \bigl[ \partial_x \tau_{12} + \partial_y \tau_{22}
+- m_y \bigl[ \partial_x \tau_{12} + \partial_y \tau_{22}
- \partial_z ( z_x \tau_{12} + z_y \tau_{22}) \bigr]
- \partial_z \tau_{23}
\label{v-disp-phys} \\
%
\partial_t W &= ...\,\,
-- m \bigl[ \partial_x \tau_{13} + \partial_y \tau_{23}
+- m_y \bigl[ \partial_x \tau_{13} + \partial_y \tau_{23}
- \partial_z ( z_x \tau_{13} + z_y \tau_{23}) \bigr]
- \partial_z \tau_{33}.
\label{w-disp-phys}
@@ -171,10 +249,10 @@
\tau_{12} &= -\mu_d K_h D_{12}
</font>
<font color="red">otag \\
%
-\tau_{13} &= -\mu_d K_h D_{13}
+\tau_{13} &= -\mu_d K_v D_{13}
</font>
<font color="red">otag \\
%
-\tau_{23} &= -\mu_d K_h D_{23}.
+\tau_{23} &= -\mu_d K_v D_{23}.
</font>
<font color="gray">otag
%
\end{align}
@@ -186,31 +264,31 @@
$D$. The continuous deformation tensor is defined as
%
\begin{align}
-D_{11} &= 2 \, m^2 \bigl[ \partial_x (m^{-1}u) - z_x \partial_z
-(m^{-1}u) \bigl]
+D_{11} &= 2 \, m_x m_y \bigl[ \partial_x (m_y^{-1}u) - z_x \partial_z
+(m_y^{-1}u) \bigl]
</font>
<font color="red">otag \\
%
-D_{22} &= 2 \, m^2 \bigl[ \partial_y (m^{-1}v) - z_y \partial_z
-(m^{-1}v) \bigl]
+D_{22} &= 2 \, m_x m_y \bigl[ \partial_y (m_x^{-1}v) - z_y \partial_z
+(m_x^{-1}v) \bigl]
</font>
<font color="black">otag \\
%
D_{33} &= 2 \, \partial_z w
</font>
<font color="red">otag \\
%
-D_{12} &= m^2 \bigl[
-\partial_y (m^{-1}u) - z_y \partial_z(m^{-1}u)
-+ \partial_x (m^{-1}v) - z_x \partial_z (m^{-1}v)
+D_{12} &= m_x m_y \bigl[
+\partial_y (m_y^{-1}u) - z_y \partial_z(m_y^{-1}u)
++ \partial_x (m_x^{-1}v) - z_x \partial_z (m_x^{-1}v)
\bigl]
</font>
<font color="red">otag \\
%
-D_{13} &= m^2 \bigl[
-\partial_x (m^{-1}w) - z_x \partial_z(m^{-1}w)
+D_{13} &= m_x m_y \bigl[
+\partial_x (m_y^{-1}w) - z_x \partial_z(m_y^{-1}w)
\bigl]
+ \partial_z (u)
</font>
<font color="red">otag \\
%
-D_{23} &= m^2 \bigl[
-\partial_y (m^{-1}w) - z_y \partial_z(m^{-1}w)
+D_{23} &= m_x m_y \bigl[
+\partial_y (m_y^{-1}w) - z_y \partial_z(m_y^{-1}w)
\bigl]
+ \partial_z (v).
</font>
<font color="gray">otag
@@ -225,12 +303,12 @@
%
\begin{align}
\partial_t (\mu_d q) = ...\,\, + \bigl[
-& \, m \bigl(\partial_x - \partial_z z_x \bigr )
-\bigl( \mu_d m K (\partial_x - z_x \partial_z ) \bigr)
+& \, m_x \bigl(\partial_x - \partial_z z_x \bigr )
+\bigl( \mu_d m_x K (\partial_x - z_x \partial_z ) \bigr)
+
</font>
<font color="gray">otag \\
-& \, m \bigl(\partial_y - \partial_z z_y \bigr )
-\bigl( \mu_d m K (\partial_y - z_y \partial_z ) \bigr)
+& \, m_y \bigl(\partial_y - \partial_z z_y \bigr )
+\bigl( \mu_d m_y K (\partial_y - z_y \partial_z ) \bigr)
+ \partial_z \mu_d K \partial_z
\bigr] q.
\label{scalar-disp-phys}
@@ -245,7 +323,7 @@
%
\begin{align}
\partial_t U &= ...\,\,
-- {\overline m}^x \bigl[ \delta_x \tau_{11} + \delta_y \tau_{12}
+- {\overline{m_x}}^x \bigl[ \delta_x \tau_{11} + \delta_y \tau_{12}
- \delta_z
( z_x {\overline{\tau_{11}}}^{x\eta} + {\overline {z_y}}^{xy}
{\overline{\tau_{12}}}^{y\eta}) \bigr]
@@ -253,7 +331,7 @@
</font>
<font color="gray">otag \\
%
\partial_t V &= ...\,\,
-- {\overline m}^y \bigl[ \delta_y \tau_{22} + \delta_x \tau_{12}
+- {\overline{m_y}}^y \bigl[ \delta_y \tau_{22} + \delta_x \tau_{12}
- \delta_z
( z_y {\overline{\tau_{22}}}^{y\eta} + {\overline {z_x}}^{xy}
{\overline{\tau_{12}}}^{x\eta}) \bigr]
@@ -261,7 +339,7 @@
</font>
<font color="gray">otag \\
%
\partial_t W &= ...\,\,
-- m \bigl[ \delta_x \tau_{13} + \delta_y \tau_{23}
+- m_x \bigl[ \delta_x \tau_{13} + \delta_y \tau_{23}
- \delta_z (
{\overline {z_x}}^{x\eta} {\overline{\tau_{13}}}^{x\eta}
+{\overline {z_y}}^{y\eta} {\overline{\tau_{23}}}^{y\eta}
@@ -290,11 +368,11 @@
</font>
<font color="red">otag \\
%
\tau_{13} &= -\overline {\mu_d}^{x}
-\overline {K_h}^{x\eta} D_{13}
+\overline {K_v}^{x\eta} D_{13}
</font>
<font color="red">otag \\
%
\tau_{23} &= -\overline {\mu_d}^{y}
-\overline {K_h}^{y\eta} D_{23},
+\overline {K_v}^{y\eta} D_{23},
</font>
<font color="gray">otag
%
\end{align}
@@ -303,37 +381,37 @@
and
%
\begin{align}
-D_{11} &= 2 \, m^2 \bigl[ \delta_x ({\overline {m^{-1}}}^x u)
-- {\overline{z_x}}^{x\eta} \delta_z ({\overline {m^{-1}}}^x u) \bigl]
+D_{11} &= 2 \, m_x m_y \bigl[ \delta_x ({\overline {m_y^{-1}}}^x u)
+- {\overline{z_x}}^{x\eta} \delta_z ({\overline {m_y^{-1}}}^x u) \bigl]
</font>
<font color="red">otag \\
%
-D_{22} &= 2 \, m^2 \bigl[ \delta_y ({\overline {m^{-1}}}^y v)
-- {\overline{z_y}}^{x\eta} \delta_z ({\overline {m^{-1}}}^y v) \bigl]
+D_{22} &= 2 \, m_x m_y \bigl[ \delta_y ({\overline {m_x^{-1}}}^y v)
+- {\overline{z_y}}^{x\eta} \delta_z ({\overline {m_x^{-1}}}^y v) \bigl]
</font>
<font color="black">otag \\
%
D_{33} &= 2 \, \delta_z w
</font>
<font color="red">otag \\
%
-D_{12} &= (\overline m^{xy})^2 \biggl[
-\delta_y ({\overline {m^{-1}}}^x u)
+D_{12} &= (\overline{m_x m_y}^{xy}) \biggl[
+\delta_y ({\overline {m_y^{-1}}}^x u)
- {\overline{z_y}}^{x\eta}
-\delta_z {\overline{({\overline {m^{-1}}}^x u)}^{y\eta}}
-+ \delta_x ({\overline {m^{-1}}}^y v)
+\delta_z {\overline{({\overline {m_y^{-1}}}^x u)}^{y\eta}}
++ \delta_x ({\overline {m_x^{-1}}}^y v)
- {\overline{z_x}}^{y\eta}
-\delta_z {\overline{({\overline {m^{-1}}}^y v)}^{x\eta}}
+\delta_z {\overline{({\overline {m_x^{-1}}}^y v)}^{x\eta}}
\biggl]
</font>
<font color="red">otag \\
%
-D_{13} &= m^2 \biggl[
-\delta_x (m^{-1} w)
+D_{13} &= m_x m_y \biggl[
+\delta_x (m_y^{-1} w)
- z_x \delta_z \overline{(m^{-1} w)}^{x\eta}
\biggl]
+ \delta_z u
</font>
<font color="red">otag \\
%
-D_{23} &= m^2 \biggl[
-\delta_y (m^{-1} w)
-- z_y \delta_z \overline{(m^{-1} w)}^{y\eta}
+D_{23} &= m_x m_y \biggl[
+\delta_y (m_y^{-1} w)
+- z_y \delta_z \overline{(m_y^{-1} w)}^{y\eta}
\biggl]
+ \delta_z v.
</font>
<font color="gray">otag
@@ -346,13 +424,13 @@
%
\begin{align}
\partial_t (\mu_d q) = ...\,\, &
-+ m \bigl[
++ m_x \bigl[
\delta_x \bigl({\overline \mu_d}^x H_1( q) \bigr)
- \mu_d \delta_z
\bigl(\overline{z_x}^x {\overline{H_1( q)}}^{x\eta} \bigr)
\bigr]
</font>
<font color="gray">otag \\
-& + m \bigl[
+& + m_y \bigl[
\delta_y \bigl({\overline \mu_d}^y H_2( q) \bigr)
- \mu_d \delta_z
\bigl(\overline{z_y}^y {\overline{H_2( q)}}^{y\eta} \bigr)
@@ -365,10 +443,10 @@
</font>
<font color="red">oindent
where
\begin{align}
-H_1( q) &= \overline{m}^x {\overline{K_h}}^x \bigl(
+H_1( q) &= \overline{m_x}^x {\overline{K_h}}^x \bigl(
\delta_x q - z_x \delta_z( {\overline{ q}}^{x\eta} )\bigr),
</font>
<font color="red">otag \\
-H_2( q) &= \overline{m}^y {\overline{K_h}}^y \bigl(
+H_2( q) &= \overline{m_y}^y {\overline{K_h}}^y \bigl(
\delta_y q - z_y \delta_z( {\overline{ q}}^{y\eta} )\bigr).
</font>
<font color="gray">otag
\end{align}
@@ -437,18 +515,14 @@
and $N$ is the Brunt-V\"ais\"al\"a frequency; the computation of $N$,
including moisture effects, is outlined in Section \ref{tke_section}.
-If $\Delta x$ is less than the user-specified critical length
-scale $l_{cr}$, then the length scale used for calculating
-both $K_h$ and $K_v$ in \eqref{k_calc} is $l_{h,v} =
-(\Delta x \Delta y \Delta z)^{1/3}$ (and $K_h = K_v = K$).
-If $\Delta x$ is greater
-than an critical length scale $l_{cr}$, then the
-length scale $l_h = \sqrt{\Delta x \Delta y}$ in the calculation of the
-horizontal eddy viscosity $K_h$ using \eqref{k_calc}, and
-$l_v = \Delta z$ for the calculation of the
+Two options are available for calculating the mixing length $l_{h,v}$ in
+\eqref{k_calc}. An isotropic length scale can be chosen where $l_{h,v}
+= (\Delta x \Delta y \Delta z)^{1/3}$ and thus $K_h = K_v = K$. The
+anisotropic option sets the horizontal mixing length $l_h = \sqrt{\Delta
+x \Delta y}$ in the calculation of the horizontal eddy viscosity $K_h$
+using \eqref{k_calc}, and $l_v = \Delta z$ for the calculation of the
vertical eddy viscosity $K_v$ using \eqref{k_calc}.
-
Additionally, the eddy viscosities
for scalar mixing are divided by the turbulent Prandtl number $P_r = 1/3$.
@@ -468,8 +542,7 @@
in this scheme), $C_k$ is a constant (typically $0.15 < C_k < 0.25$),
and $l$ is a length scale.
-If the grid spacing $\Delta x$ is less than the critical length scale
-$l_{cr}$, then
+If the isotropic mixing option is chosen,
%
\begin{align}
l_{h,v} & = \hbox{min}
@@ -492,8 +565,8 @@
horizontal and vertical eddy viscosities
are equivalent.
-If the grid spacing $\Delta x$ is greater than the critical length scale
-$l_{cr}$, then $l_h=\sqrt{\Delta x \Delta y}$ for the calculation of $K_h$.
+If the anisotropic mixing option is chosen,
+then $l_h=\sqrt{\Delta x \Delta y}$ for the calculation of $K_h$.
For calculating $K_v$,
%
\begin{align} l_v &= \hbox{min}
@@ -597,8 +670,7 @@
\subsubsection{Dissipation}
-If $\Delta x$ is less than the critical length scale $l_{cr}$, the
-dissipation term in \eqref{tke} is
+The dissipation term in \eqref{tke} is
%
\begin{equation}
\hbox{dissipation} = - {Ce^{3/2} \over l},
@@ -622,34 +694,36 @@
</font>
<font color="red">otag
\end{equation}
-%%%%%%%%%%%%%%%
+\subsection{Sixth-Order Spatial Filter on Coordinate Surfaces}
-If $\Delta x$ is greater than the critical length scale $l_{cr}$, the
-dissipation term in \eqref{tke} is
+A sixth order spatial filter is available that is applied on horizontal
+coordinate surfaces.
+The diffusion scheme is that proposed by \cite{xue.2000}. Its application
+in the ARW is described by \cite{knievel.et.al.2007}. The filter can be
+expressed as
%
\begin{equation}
-\hbox{dissipation} = - { 2 \sqrt{2} \over 15}
-{e^{3/2} \over l},
-</font>
<font color="blue">otag
+\partial_t (\mu_d a) = ... {\beta \, 2^{-6} \over 2 \Delta t}
+\left[ \Delta x^6 \delta_x (\overline{\mu_d}^x F_x) + \Delta y^6 \delta_y (\overline{\mu_d}^y F_y)\right]
+\label{sixth_order_diffusion}
\end{equation}
%
</font>
<font color="red">oindent
-where
-%
-\begin{equation}
-l = {kz \over 1 + kz/l_0},
-</font>
<font color="red">otag
-\end{equation}
-%
-\begin{equation}
-l_0 = min\biggl( {\alpha_b \int_0^{z_i} \sqrt{e}z\,dz
-\over \int_0^{z_i} \sqrt{e}\,dz}, \,\, 80 \biggr),
-</font>
<font color="red">otag
-\end{equation}
-%
-</font>
<font color="gray">oindent
-$\alpha_b = 0.2$, and $k = 0.4$ is the von Karman constant.
+The diffusive fluxes $F_x = \delta_x^5 (a)$ and $F_y = \delta_y^5 (a)$.
+$\beta$ is a user-specified filtering coefficient representing the
+damping applied to $2\Delta$ waves for each filter application.
+The user can choose to enforce monotonicity in the filtering. In the
+monotonic option, the diffusive fluxes are set to zero if up-gradient
+diffusion is indicated -- $F_x = 0$ if $F_x \, \delta_x(a) < 0$ and $F_y
+= 0$ if $F_y \, \delta_y(a) < 0$. Note that the filter is applied in
+computational space; map factors are not taken into account
+and the filter parameter $\beta$ is dimensionless. Thus, the filter
+scales with the timestep and gridsize and it does not conserve scalar
+mass. This explicit diffusion acts on all three components of wind,
+potential temperature, all moisture variables and passive scalars, and
+subgrid turbulence kinetic energy. Further details can be found in
+\cite{knievel.et.al.2007}.
%%%%%%%%%%%%%%%%%%%%%%%%%%
%
@@ -657,7 +731,6 @@
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
-
\section{Filters for the Time-split RK3 scheme}
Three filters are used in the ARW time-split RK3
@@ -668,7 +741,7 @@
equation and geopoential equation. Each of these is described in the
following sections.
-\subsection{Divergence Damping}
+\subsection{Three-Dimensional Divergence Damping}
The damping of the full mass divergence is a filter for acoustic
modes in the ARW solver. This 3D mass divergence damping is described
@@ -736,7 +809,7 @@
Forward-in-time weighting of the vertically-implicit
acoustic-time-step terms damps instabilities associated
vertically-propagating sound waves. The forward weighting also damps
-instabilities associated with sloping mode levels and horizontally
+instabilities associated with sloping model levels and horizontally
propagating sound waves
\citep[see][]{durran_klemp83, dudhia95}. The off-centering
is accomplished by using a positive (non-zero) coefficient $\beta$
@@ -746,18 +819,22 @@
typically used in the ARW applications, independent of time step or
grid size.
-\section{Other Damping}
+\section{Formulations for Gravity Wave Absorbing Layers}
-\subsection{Gravity Wave Absorbing Layer}
+There are three formulations for absorbing vertically-propagating
+gravity waves so as to prevent unphysical wave reflection off the
+domain upper boundary from contaminating ARW solutions.
-A gravity-wave absorbing layer is available in the ARW solver.
-The absorbing layer increases the second-order horizontal and
-vertical eddy viscosities in the absorbing layer using the
-following formulation:
+\subsection{Absorbing Layer Using Spatial Filtering}
+\label{absorbing_layer_spatial}
+
+This formulation of the absorbing layer increases the second-order
+horizontal and vertical eddy viscosities in the absorbing layer using
+the following formulation:
%
\begin{equation}
K_{dh} = {\Delta x^2\over \Delta t} \gamma_g
-\cos \biggl({z_{top} - z \over z_{d}}{\pi \over 2} \biggr),
+\cos \biggl({\pi \over 2} \,{z_{top} - z \over z_{d}} \biggr),
</font>
<font color="gray">otag
\end{equation}
%
@@ -766,7 +843,7 @@
%
\begin{equation}
K_{dv} = {\Delta z^2\over \Delta t} \gamma_g
-\cos \biggl({z_{top} - z \over z_{d}}{\pi \over 2} \biggr).
+\cos \biggl({\pi \over 2} \, {z_{top} - z \over z_{d}} \biggr).
</font>
<font color="gray">otag
\end{equation}
%
@@ -782,36 +859,91 @@
maximum of ($K_h$, $K_{dh}$) and ($K_v$, $K_{dv}$).
The effect of this filter on gravity waves is discussed in
\citet{klemp_and_lilly78}, where guidance on the choice of
-the damping coefficeint $\gamma_g$ can also be found.
+the damping coefficient $\gamma_g$ can also be found.
-\subsection{Rayleigh Damping Layer}
+\subsection{Implicit Rayleigh Damping for the Vertical Velocity}
+\label{rayleigh-w}
-A Rayleigh damping layer is also available in the ARW solver. This
+Implicitly damping the vertical velocity within the implicit solution
+algorithm for the vertically propagating acoustic modes has been found
+to be a very effective and robust absorbing layer formulation by
+\citet{klemp_et_al_2008}. We recommend its use in both large-scale and
+small-scale applications, and in idealized and real-data applications.
+It has proven more effective than the spatial-filter-based absorbing
+layer described in Section \ref{absorbing_layer_spatial}, and it is more
+effective and more generally applicable than traditional Rayleigh
+damping because of its implicit nature and the fact that it does not
+need a reference state. This formulation has been introduced in WRFV3.
+It can only be used with the nonhydrostatic dynamics option.
+
+In the vertically-implicit solution procedure for $W''^{\tau+\Delta
+\tau}$ and $\phi''^{\tau+\Delta \tau}$ in the acoustic step (step 5 in
+Figure \ref{time_integration_figure}), equation \eqref{geo-small-step}
+is used to eliminate $\phi''^{\tau+\Delta \tau}$ from
+\eqref{w-small-step} after which the tridiagonal equation in the
+vertical direction for
+$W''^{\tau+\Delta \tau}$ is solved. Afterwards $\phi''^{\tau+\Delta \tau}$ is
+recovered using \eqref{geo-small-step}. In the solution procedure that
+includes the implicit Rayleigh damping for $W$, after the tridiagonal
+equation for $W''$ %$W''^{\tau+\Delta \tau}$
+
+is solved and before the
+geopotential $\phi$ is updated, the implicit Rayleigh damping is
+included by calculating $W''^{\tau+\Delta \tau}$ using
+%
+\begin{equation}
+W''^{\tau+\Delta \tau} = \tilde{W}''^{\tau+\Delta \tau}
+-\tau(z) \Delta \tau W''^{\tau+\Delta \tau}
+\label{w-rayleigh}
+\end{equation}
+%
+where $\tilde{W}''^{\tau+\Delta \tau}$ is the solution to
+\eqref{w-small-step}. The geopotential is then updated as usual using
+\eqref{geo-small-step} with the updated damped vertical velocity from
+\eqref{w-rayleigh}.
+The perturbation pressure and specific volume
+are computed as before using \eqref{p-linear} and \eqref{small-hydro}.
+
+The variable $\tau(z)$ defines the vertical structure of the damping
+layer, and has a form similar to the Rayleigh damper in
+\citet{durran_klemp83} and also used in the traditional Rayleigh damping
+formulation discussed in Section \ref{traditional-rayleigh}
+%
+\[ \tau (z) = \left\{ \begin{array}{cc}
+ \gamma_r\sin^2 \left[ \frac{\pi}{2}
+\left( 1 - \frac{z_{top}-z}{z_d} \right) \right] &
+\mbox{for $z \geq (z_{top}-z_d)$};\\
+ 0 & \mbox{otherwise}, \end{array} \right. \]
+%
+where $\gamma_r$ is a user specified damping coeficient, $z_{top}$ is
+the height of the model top for a particular grid column, and $z_d$ is
+the depth of the damping layer (from the model top). A typical value
+for the damping coefficient for this formulation is
+$\gamma_r = 0.2 \,s^{-1}$ ($\sim 10N$ for typical stratospheric values of
+the Brunt-V\"ais\"all\"a frequency).
+A complete
+analysis of this filter can be found in \citet{klemp_et_al_2008}.
+
+\subsection{Traditional Rayleigh Damping Layer}
+\label{traditional-rayleigh}
+
+A traditional Rayleigh damping layer is also available in the ARW solver. This
scheme applies a tendency to $u$, $v$, $w$, and $\theta$ to gradually
relax the variable back to a predetermined reference state value,
%
\begin{eqnarray*}
- \frac{\partial u }{\partial t} & = & \tau (z) \left( u - \overline{u} \right) ,\\
- \frac{\partial v }{\partial t} & = & \tau (z) \left( v - \overline{v} \right) ,\\
- \frac{\partial w }{\partial t} & = & \tau (z) w ,\\
- \frac{\partial \theta}{\partial t} & = & \tau (z) \left( \theta - \overline{\theta} \right).
+ \frac{\partial u }{\partial t} & = & -\tau (z) \left( u - \overline{u} \right) ,\\
+ \frac{\partial v }{\partial t} & = & -\tau (z) \left( v - \overline{v} \right) ,\\
+ \frac{\partial w }{\partial t} & = & -\tau (z) w ,\\
+ \frac{\partial \theta}{\partial t} & = & -\tau (z) \left( \theta - \overline{\theta} \right).
\end{eqnarray*}
%
Overbars indicate the reference state fields, which are functions of $z$
only and are typically defined as the initial horizontally homogeneous
fields in idealized simulations. The reference state vertical velocity
-is assumed to be zero. The variable $\tau$ defines the vertical
-structure of the damping layer, and has a form similar to the Rayleigh
-damper in \citet{durran_klemp83},
-\[ \tau (z) = \left\{ \begin{array}{cc}
- -\gamma_r\sin^2 \left[ \frac{\pi}{2}
-\left( 1 - \frac{z_{top}-z}{z_d} \right) \right] &
-\mbox{for $z \geq (z_{top}-z_d)$};\\
- 0 & \mbox{otherwise}, \end{array} \right. \]
-where $\gamma_r$ is a user specified damping coeficient, $z_{top}$ is
-the height of the model top for a particular grid column, and $z_d$ is
-the depth of the damping layer (from the model top).
-
+is assumed to be zero. The vertical structure of the damping layer is
+the same as that used for the implicit vertical velocity damping
+described in Section \ref{rayleigh-w}.
Because the model surfaces change height with time in the ARW solver, the
reference state values at each grid point need to be recalculated at
every time step. Thus, a linear interpolation scheme is used to
@@ -822,6 +954,8 @@
\citet{klemp_and_lilly78}, where guidance on the choice of
the damping coefficeint $\gamma_r$ can also be found.
+\section{Other Damping}
+
\subsection{Vertical-Velocity Damping}
This is also called $w$-damping.
</font>
</pre>