Vertical discretization

The model uses the so-called Lorenz vertically staggered grid which is shown in Fig. 8 for a $N$ layer model ($N+1$ levels): the index is $k=1$ at the top of the atmosphere, and increases towards the surface (with increasing pressure). Mass and wind component variables are defined at layer centers, while the vertical velocities are defined at the levels or layer interfaces.

Figure 8: The model vertical discretization. A hydrostatic pressure-based vertical coordinate is assumed. Atmospheric state variables are defined at layer centers (dashed lines), while the vertical advection and turbulent fluxes, $F$, are defined along layer interfaces (solid lines: pressure levels). ${\tilde \phi }$ indicates a layer variable $\phi$ averaged to a level, and ${\overline \vartheta }$ indicates a level variable $\vartheta$ interpolated to a layer center. The bottom boundary of the grid is along the terrain-following surface hydrostatic pressure, $\pi _s$, which corresponds to the geopotential height, $\Phi _s$ (or equivalently, $\Phi _{N+1}$).

The pressure along each hybrid ($\eta$) surface is simply

$$\pi_{d,k} = A_k + B_k \pi_{d,s} \hskip.5in \left( k=1,N+1\right)$$ (159)

The specific density (vertical pressure derivative) is simply discretized as

$$\mu_{d,k} = {\frac{\Delta A_{\eta,\,k} + \Delta B_{\eta,\,k} \, \pi_{d,s} }{ \Delta\eta_{k} }}$$ (160)

where the vertical pressure and coordinate differences are

$$\Delta\eta_k = \eta_{k+1}-\eta_{k}$$ (161)

$$\Delta \pi_{d,k} = \pi_{d,k+1}-\pi_{d,k}$$ (162)

The finite difference form of the surface pressure tendency (Eq. 34), which is simply the sum of the column mass convergence, can then be readily obtained:

$${\frac{\partial \pi_{d,s} }{\partial t}} = - \sum_{k=1}^{N} \nabla\cdot{\left( {\bf V}_k\mu_{d,k}\right)} \Delta\eta_{k}$$ (163)

The vertical mass coordinate velocity (Eq.s 51-52) can then be discretized as

$${\dot\pi}_k = {\dot\eta}_k \mu_{d,k} = - B_{\eta,\,k} {\frac{\par... ...}{\partial t}} - \sum_{l=1}^k \nabla \cdot ({\bf V}_l \mu_{d,l}) \,\Delta\eta_l$$ (164)

subject to the boundary conditions

$${\dot\pi}_1 = {\dot\pi}_{N+1} = 0$$ (165)

The discretized vertical advection terms for the heat, the $u$ and $v$ wind components and the remaining scalars are then defined as, respectively:

$${\frac{\partial {\left( {\dot\eta} \mu \theta \right)}_k }{\partial\eta}} = {\frac{1 }{\Delta\eta_k}} \left( {\dot\pi}_{k+1} {\hat{\theta}}_{k+1} - {\dot\pi}_{k} {\hat{\theta}}_{k} \right)$$ (166)

$${\frac{\partial {\left( {\dot\eta} \mu u \right)}_k }{\partial\eta}} = {\frac{1 }{\Delta\eta_k}} \left( {\overline{\dot\pi}^u}_{k+1} {\hat{u}}_{k+1} - {\overline{\dot\pi}^u}_{k} {\hat{u}}_{k} \right)$$ (167)

$${\frac{\partial {\left( {\dot\eta} \mu v \right)}_k }{\partial\eta}} = {\frac{1 }{\Delta\eta_k}} \left( {\overline{\dot\pi}^v}_{k+1} {\hat{v}}_{k+1} - {\overline{\dot\pi}^v}_{k} {\hat{v}}_{k} \right)$$ (168)

$${\frac{\partial {\left( {\dot\eta} \mu w \right)}_k }{\partial\eta}} = {\frac{1 }{\Delta\eta_k}} \left( {\overline{\dot\pi}}_{k+1} {\overline{w}}_{k+1} - {\overline{\dot\pi}}_{k} {\overline{w}}_{k} \right)$$ (169)

$${\frac{\partial {\left( {\dot\eta} \mu \varphi \right)}_k }{\partial\eta}} = {\frac{1 }{\Delta\eta_k}} \left( {\dot\pi}_{k+1} {\hat\varphi}_{k+1} - {\dot\pi}_{k} {\hat\varphi}_{k} \right)$$ (170)

where ${\overline{\phi}}^u$ and ${\overline{\phi}}^v$ represent the horizontal interpolation of variable $\phi$ from the horizontal mass center of the grid cell to the $u$ and $v$ points on the staggered C-grid, respectively, and ${\overline{\phi}}$ represents the vertical interpolation from (hydrostatic pressure) levels to the vertical mass center of the grid cell. The level values of the aforementioned quantities are interpolated from the layer centers to levels. For variables defined on levels which must be interpolated to layer centers, the layer average value is simply the average of the surrounding level values:

$${\overline{\phi}}_{k} = \frac{\phi_k + \phi_{k+1}}{2}$$ (171)

which is used for $w$ and ${\overline{\dot\pi}}^\eta$. The potential temperature is interpolated using

$${\hat\theta}_{k} = \frac{ \left( \Pi_{\pi,d_k} - {\overline{\Pi}}... ...heta_{k} }{ {\overline{\Pi}}_{\pi d_{k}} \,-\, {\overline{\Pi}}_{\pi d_{k-1}} }$$ (172)

where ${\Pi_\pi}$ represents the Exner function along the sigma surfaces, which is computed using the hydrostatic pressure, $\pi$, from Eq. 15 as

$$\Pi_{\pi d_{k}} = {\left({\frac{\pi_{d k}}{p_0}}\right)}^{R/C_{pd}} \hskip1.in (k=1,N+1)$$ (173)

The Exner function at the layer centers (corresponding to the location of the prognostic variables) is defined as

$${\overline{\Pi}}_{\pi d_{k}} = \frac{1}{1+\left(R/C_{pd}\right)} ... ...pi_{d k} \Pi_{\pi d_{k}}} {\pi_{d k+1} - \pi_{d k} } \right] \hskip1.in (k=1,N)$$ (174)

The consistent layer average hydrostatic pressure, ${\overline\pi}_k$, can computed from Eq. 174 but using ${\overline{\Pi}}_{\pi_{k}}^\eta$.

Because of the rapid decrease of saturation vapor pressure with increasing height and the relatively coarse vertical layering (especially higher up in the atmosphere), $q$ is interpolated to levels using

$${\hat q}_k = {\frac {\varsigma_k + \varsigma_{k-1}} {\left(\varsigma_k/q_k\right) + \left(\varsigma_{k-1}/q_{k-1}\right)} }$$ (175)

where we have chosen $\varsigma_k=\Delta\Pi_{\pi_k}$. This method was found to result in good relative humidity values, especially at the upper levels in the atmosphere, but it should be noted that numerous other methods exist. The key is to account for the potential rapid decrease in $q$ (owing to the strong temperature and pressure dependence of $q_{sat}$), especially for relatively coarse (vertical) grids.

The interpolant for the remaining variables denoted by $\phi$ is simply:

$${\hat{\phi}}_k = \frac{ \Delta\eta_{k-1}\phi_{k} \,+\, \Delta\eta... ..._{k-1}\,+\, \Delta\eta_{k}} \hskip.5in \left( \phi = \, u,\, v,\,\varphi\right)$$ (176)

where $\varphi$ represents any additional scalars (cloud water, TKE, etc.).

The vertical coordinate velocities are interpolated from mass to velocity points using the usual operators as

$${\overline{\dot\pi}^x}_{i,j,k} = \left({\dot\pi}_{i,j,k} + {\dot\pi}_{i-1,j,k}\right)/2$$ (177)

$${\overline{\dot\pi}^y}_{i,j,k} = \left({\dot\pi}_{i,j,k} + {\dot\pi}_{i,j-1,k}\right)/2$$ (178)

The non-hydrostatic surface pressure (Eq. 56) can then simply be expressed as

$$p_s ={\hat{p}}_{N+1} \,=\, p_1 \,+\, \sum_{k=1}^N \left({\overline{\epsilon}}_{k} + 1 \right)\, \frac{\mu_{d k}}{\varphi_{q k}} \Delta\eta_k$$ (179)

where $\epsilon$ is vertically interpolated to layer centers using Eq. 172. Note that $p_1=\pi_{d 1}$. In Eq. 180, the layer centered value of $\epsilon$ is simply the average of the interfacial values above and below. This method ignores the fact that the non-hydrostatic pressure is interpolated using the Exner function to layer centers, but the approximation turns out to be good (through testing).