3.2 Semi-Lagrangian Dynamical Core

3.2.1 Introduction

The two-time-level semi-implicit semi-Lagrangian spectral transform dynamical core in CAM 3.0 evolved from the three-time-level CCM2 semi-Lagrangian version detailed in Williamson and Olson [185] hereafter referred to as W&O94. As a first approximation, to convert from a three-time-level scheme to a two-time-level scheme, the time level index n-1 becomes n, the time level index n becomes n+ $\frac{1}{ 2}$, and $2 \Delta t$ becomes $\Delta t$. Terms needed at n+ $\frac{1}{ 2}$ are extrapolated in time using time n and n-1 terms, except the Coriolis term which is implicit as the average of time n and n+1. This leads to a more complex semi-implicit equation to solve. Additional changes have been made in the scheme to incorporate advances in semi-Lagrangian methods developed since W&O94. In the following, reference is made to changes from the scheme developed in W&O94. The reader is referred to that paper for additional details of the derivation of basic aspects of the semi-Lagrangian approximations. Only the details of the two-time-level approximations are provided here.

3.2.2 Vertical coordinate and hydrostatic equation

The semi-Lagrangian dynamical core adopts the same hybrid vertical coordinate ($\eta$) as the Eulerian core defined by

$$p(\eta, p_s) = A (\eta)p_o + B(\eta) p_s \; ,$$ (3.263)

where $p$ is pressure, $p_s$ is surface pressure, and $p_o$ is a specified constant reference pressure. The coefficients $A$ and $B$ specify the actual coordinate used. As mentioned by Simmons and Burridge [157] and implemented by Simmons and Strüfing [158] and Simmons and Strüfing [159], the coefficients $A$ and $B$ are defined only at the discrete model levels. This has implications in the continuity equation development which follows.

In the $\eta$ system the hydrostatic equation is approximated in a general way by

$$\Phi_k = \Phi_s + R \sum^K_{l=k} H_{kl} (p) T_{vl}$$ (3.264)

where k is the vertical grid index running from 1 at the top of the model to $K$ at the first model level above the surface, $\Phi_k$ is the geopotential at level $k$, $\Phi_s$ is the surface geopotential, $T_v$ is the virtual temperature, and R is the gas constant. The matrix $H$, referred to as the hydrostatic matrix, represents the discrete approximation to the hydrostatic integral and is left unspecified for now. It depends on pressure, which varies from horizontal point to point.

3.2.3 Semi-implicit reference state

The semi-implicit equations are linearized about a reference state with constant $T^r$ and $p_s^r$. We choose

$$T^r = 350 {\rm K},      p_s^r = 10^5 {\rm Pa}$$ (3.265)

3.2.4 Perturbation surface pressure prognostic variable

To ameliorate the mountain resonance problem, Ritchie and Tanguay [148] introduce a perturbation ${\rm ln} p_{s}$ surface pressure prognostic variable

$${\rm ln} p'_{s} = {\rm ln} p_{s} - {\rm ln} p_{s}^*$$ (3.266)

$${\rm ln} p_{s}^* = - \frac{\Phi_s}{R T^r}$$ (3.267)

The perturbation surface pressure, ${\rm ln} p'_{s}$, is never actually used as a grid point variable in the CAM 3.0 code. It is only used for the semi-implicit development and solution. The total ${\rm ln} p_{s}$ is reclaimed in spectral space from the spectral coefficients of $\Phi_s$ immediately after the semi-implicit equations are solved, and transformed back to spectral space along with its derivatives. This is in part because $\nabla ^4{\rm ln} p_{s}$ is needed for the horizontal diffusion correction to pressure surfaces. However the semi-Lagrangian CAM 3.0 default is to run with no horizontal diffusion.

3.2.5 Extrapolated variables

Variables needed at time ( $n+\frac{1}{2}$) are obtained by extrapolation

$$\left(   \right)^{n+\frac{1}{2}} = \frac{3}{2} \left(   \right)^{n} - \frac{1}{2} \left(   \right)^{n-1}$$ (3.268)

3.2.6 Interpolants

Lagrangian polynomial quasi-cubic interpolation is used in the prognostic equations for the dynamical core. Monotonic Hermite quasi-cubic interpolation is used for tracers. Details are provided in the Eulerian Dynamical Core description. The trajectory calculation uses tri-linear interpolation of the wind field.

3.2.7 Continuity Equation

The discrete semi-Lagrangian, semi-implicit continuity equation is obtained from (16) of W&O94 modified to be spatially uncentered by a fraction $\epsilon$, and to predict ${\rm ln} p_{s}'$

$$\Delta B_{_{l}} \left\{ \left({\rm ln} p_{s{_{l}}}'\right)^{n+1}_{_{A}} - \left[... ...)^{n} + { \Phi_s \over RT^r } \right] _{D_{2}} \right\} \Bigg/ \Delta t =$$

$$- {1\over 2} \bigg\lbrace \left[ \left(1 + \epsilon \right) \Delt... ...dot\eta {\partial p \over \partial\eta}\right)_l\right]^{n}_{D_{2}}\bigg\rbrace$$ (3.269)

$$-\left({1 \over p_s} \delta_{_{l}} \Delta p_{_{l}} \right)^{n+{1\... ...ft( {\boldsymbol {V}}_{_{l}} \cdot \nabla \;\Phi_s \right)^{n+{1\over 2}}_{M_2}$$

$$-\bigg\lbrace{1 \over 2} \left[ \left(1 + \epsilon \right) \left(... ...r_s} \delta_{_{l}} \Delta p_{_{l}}^r \right)^{n+{1\over 2}}_{M_{2}}\bigg\rbrace$$

where

$$\Delta (    )_l = (    )_{l+ {1 \over 2}} - (    )_{l - {1 \over 2}}$$ (3.270)
and

$$\left(     \right)^{n+{1\over 2}}_{M_{2}} = {1\over 2} \left[ \left(1 + \epsilon \right)\left(     \right)^... ...} + \left(1 - \epsilon \right)\left(     \right)^{n+{1\over 2}}_{D_{2}} \right]$$ (3.271)

$\Delta (    )_l$ denotes a vertical difference, $l$ denotes the vertical level, $A$ denotes the arrival point, $D_2$ the departure point from horizontal (two-dimensional) advection, and $M_2$ the midpoint of that trajectory.

The surface pressure forecast equation is obtained by summing over all levels and is related to (18) of W&O94 but is spatially uncentered and uses ${\rm ln} p_{s}'$

$$\left({\rm ln} p'_s \right)^{n+1}_{_{A}} = \sum^K_{l=1} \Delta B_l \left[ \left({\rm ln} p_{s{_{l}}} \right... ...\over p_s} \dot\eta {\partial p \over \partial\eta}\right)_l\right]^{n}_{D_{2}}$$

$$\phantom{=} - \Delta t \sum^K_{l=1}\left({1 \over p_s} \delta_l \Delta p_{_{l... ...ft( {\boldsymbol {V}}_{_{l}} \cdot \nabla \;\Phi_s \right)^{n+{1\over 2}}_{M_2}$$ (3.272)

$$\phantom{=} - \Delta t \sum^K_{l=1} {1 \over p^r_s} \bigg\lbrace{1 \over 2} \... ... \left(\delta_l\right)^{n+{1\over 2}}_{_{M_{2}}} \bigg\rbrace \Delta p^r_{_{l}}$$

The corresponding $\left({1 \over p_s} \dot\eta {\partial p \over \partial \eta}\right)$ equation for the semi-implicit development follows and is related to (19) of W&O94, again spatially uncentered and using ${ln} p_{s}'$.

$$\left(1 + \epsilon \right) \left({1 \over p_s} \dot\eta {\partial p \over \partial \eta}\right)^{n+1}_{k+{1 \over 2}} = - {2 \over \Delta t} \bigg\lbrace B_{k+ {1\over 2}} \left({\rm ln... ...m ln}\; p_{s_l} \right)^{n} + { \Phi_s \over RT^r } \right] _{D_2} \bigg\rbrace$$

$$- \sum^k_{l=1}\left[ \left(1 - \epsilon \right)\Delta \left({1 \over p_s} \dot\eta {\partial p \over \partial\eta}\right)_l \right]^{n}_{D_{2}}$$ (3.273)

$$- 2 \sum^k_{l=1} \left({1 \over p_s} \delta_l \Delta p_l \right)^... ...ft( {\boldsymbol {V}}_{_{l}} \cdot \nabla \;\Phi_s \right)^{n+{1\over 2}}_{M_2}$$

$$- 2 \sum^k_{l=1} {1 \over p^r_s} \bigg\lbrace{1 \over 2} \left[ \... ... \left(\delta_l\right)^{n+{1\over 2}}_{_{M_{2}}} \bigg\rbrace \Delta p^r_{_{l}}$$

This is not the actual equation used to determine $\left({1 \over p_s} \dot\eta {\partial p \over \partial \eta}\right)$ in the code. The equation actually used in the code to calculate $\left({1 \over p_s} \dot\eta {\partial p \over \partial \eta}\right)$ involves only the divergence at time ($n+1$) with $\left({\rm ln}\; p'_s \right)^{n+1}$ eliminated.

$$\left(1 + \epsilon \right) \left({1 \over p_s} \dot\eta {\partial p \over \partial \eta}\right)^{n+1}_{k+{1 \over 2}} =$$

$${2 \over \Delta t} \left[ \sum^k_{l=1} -  B_{k+{1\over 2}}\sum^K_{l=1}\right] \Delt... ...ft[ \left({\rm ln}\; p_{s_l} \right)^{n} + { \Phi_s \over RT^r } \right] _{D_2}$$

$$- \left[ \sum^k_{l=1} -  B_{k+{1\over 2}}\sum^K_{l=1}\right] \left... ... \over p_s} \dot\eta {\partial p \over \partial \eta}\right)_l\right]^{n}_{D_2}$$

$$- 2 \left[ \sum^k_{l=1} -  B_{k+{1\over 2}}\sum^K_{l=1}\right] \left({1 \over p_s} \delta_l \Delta p_l \right)^{n+{1\over 2}}_{M_2}$$ (3.274)

$$+ 2 \left[ \sum^k_{l=1} -  B_{k+{1\over 2}}\sum^K_{l=1}\right] {\Del... ...ft( {\boldsymbol {V}}_{_{l}} \cdot \nabla \;\Phi_s \right)^{n+{1\over 2}}_{M_2}$$

$$- 2 \left[ \sum^k_{l=1} -  B_{k+{1\over 2}}\sum^K_{l=1}\right] {1 \o... ...ght] - \left(\delta_l \right)^{n+{1\over 2}}_{_{M_2}} \bigg\rbrace \Delta p^r_l$$

The combination $\left[ \left({\rm ln} p_{s{_{l}}} \right)^{n} + { \Phi_s \over RT^r } + {1 \o... ...{\boldsymbol {V}} \cdot \nabla \;\Phi_s \right)^{n+{1\over 2}} \right] _{D_{2}}$ is treated as a unit, and follows from (3.271).

3.2.8 Thermodynamic Equation

The thermodynamic equation is obtained from (25) of W&O94 modified to be spatially uncentered and to use ${\rm ln} p_{s}'$. In addition Hortal's modification [172] is included, in which

$${d \over dt} \left[ -\left(p_s B {\partial T \over \partial p}\right)_{ref} { \Phi_s \over RT^r } \right]$$ (3.275)

is subtracted from both sides of the temperature equation. This is akin to horizontal diffusion which includes the first order term converting horizontal derivatives from eta to pressure coordinates, with $\left({\rm ln}\;p_s\right)$ replaced by $-{ \Phi_s \over RT^r }$, and $\left(p_s B {\partial T \over \partial p}\right)_{ref}$ taken as a global average so it is invariant with time and can commute with the differential operators.

$${T_A^{n+1} - T_D^{n} \over \Delta t} = \left\{ \left\{\left[ -\left(p_s B(\eta) {\partial T \over \parti... ...ght)_{ref} { \Phi_s \over RT^r } \right]_D^{n} \right\} \Bigg/ \Delta t \right.$$

$$\phantom{=} \left. +{ 1 \over RT^r }\left[ \left(p_s B(\eta) {\partial T \ove... ...a) {\partial T \over \partial p}\right)_{ref} \right]_M^{n+{1\over 2}} \right\}$$

$$\phantom{=} +\left({RT_v \over c_p^*} {\omega \over p} \right)^{n+{1\over 2}}_M + Q^{n}_M$$

$$\phantom{=} + {RT^r \over c_p} {p_s^r \over p^r}\left[ B(\eta) {d_2\;{\rm ln}... ...eft({1 \over p_s} \dot \eta {\partial p \over \partial \eta} \right)}^t \right]$$ (3.276)

$$\phantom{=} - {RT^r \over c_p} {p_s^r \over p^r} \left[\left({p \over p_s}\right) \left( {\omega \over p} \right) \right]^{n+{1\over 2}}_M$$

$$\phantom{=} - {RT^r \over c_p} {p_s^r \over p^r} B(\eta) \left[ {1 \over RT^r} {\boldsymbol {V}} \cdot \nabla \;\Phi_s \right]^{n+{1\over 2}}_{M_2}$$

Note that $Q^{n}$ represents the heating calculated to advance from time $n$ to time $n+1$ and is valid over the interval.

The calculation of $\left(p_s B {\partial T \over \partial p}\right)_{ref}$ follows that of the ECMWF (Research Manual 3, ECMWF Forecast Model, Adiabatic Part, ECMWF Research Department, 2nd edition, 1/88, pp 2.25-2.26) Consider a constant lapse rate atmosphere

$$T = T_0\left({p\over p_0}\right)^{R \gamma / g}$$ (3.277)

$${\partial T \over \partial p} = {1 \over p}{R \gamma \over g} T_0\left({p\over p_0}\right)^{R \gamma / g}$$ (3.278)

$$p_s B {\partial T \over \partial p} = B {p_s \over p}{R \gamma \over g} T$$ (3.279)

$$\left( p_s B {\partial T \over \partial p} \right)_{ref} = B_k {(p_s)_{ref} \over (p_k)_{ref}} {R \gamma \over g} (T_k)_{ref}  \hbox{for}   (T_k)_{ref} > T_C$$ (3.280)

$$\left( p_s B {\partial T \over \partial p} \right)_{ref} = 0  for   (T_k)_{ref} \leq T_C$$ (3.281)

$$(p_k)_{ref} = A_k p_0 + B_k (p_s)_{ref}$$ (3.282)

$$(T_k)_{ref} = T_0\left({(p_k)_{ref} \over (p_s)_{ref} }\right)^{R \gamma / g}$$ (3.283)

$$(p_s)_{ref} = 1013.25 {\rm mb}$$ (3.284)

$$T_0 = 288 {\rm K}$$ (3.285)

$$p_0 = 1000 {\rm mb}$$ (3.286)

$$\gamma = 6.5 {\rm K / km}$$ (3.287)

$$T_C = 216.5 {\rm K}$$ (3.288)

3.2.9 Momentum equations

The momentum equations follow from (3) of W&O94 modified to be spatially uncentered, to use ${\rm ln} p_{s}'$, and with the Coriolis term implicit following Côté and Staniforth [43] and Temperton [171]. The semi-implicit, semi-Lagrangian momentum equation at level $k$ (but with the level subscript $k$ suppressed) is

$${{\boldsymbol {V}}^{n+1}_{_{A}} - {\boldsymbol {V}}^{n}_{_{D}} \over \Delta t} = -{1 \over 2}\bigg\lbrace \left(1 + \epsilon \right) \left[f{\b... ...\times {\boldsymbol {V}} \right]^{n}_{_{D}}\bigg\rbrace + {\boldsymbol {F}}^n_M$$

$$\phantom{=} -{1 \over 2}\bigg\lbrace \left(1 + \epsilon \right) \left[ \nab... ...ght) + RT_v {B \over p} p_s \nabla {\rm ln} p_s \right]^{n+{1\over 2}}_{_{A}}$$

$$\phantom{=} \phantom{-{1 \over 2}} + \left(1 - \epsilon \right) \left[ \nabla... ...{B \over p} p_s \nabla {\rm ln} p_s \right]^{n+{1\over 2}}_{_{D}}\bigg\rbrace$$

$$\phantom{=} - {1 \over 2} \bigg\lbrace \left(1 + \epsilon \right) \nabla \l... ...ol {H}}_k^r \cdot {\boldsymbol {T}} + RT^r {\rm ln} p_s' \right]^{n+1}_{_{A}}$$ (3.289)

$$\phantom{=} \phantom{- {1 \over 2}} - \left(1 + \epsilon \right) \nabla \left... ..._k^r \cdot {\boldsymbol {T}} + RT^r {\rm ln} p_s\right]^{n+{1\over 2}}_{_{A}}$$

$$\phantom{=} \phantom{- {1 \over 2}} + \left(1 - \epsilon \right) \nabla\left[... ...symbol {H}}_k^r \cdot {\boldsymbol {T}} + RT^r {\rm ln} p_s\right]^{n}_{_{D}}$$

$$\phantom{=} \phantom{- {1 \over 2}}- \left(1 - \epsilon \right) \nabla\left[\... ...\boldsymbol {T}} + RT^r {\rm ln} p_s\right]^{n+{1\over 2}}_{_{D}}\bigg\rbrace$$

The gradient of the geopotential is more complex than in the $\sigma$ system because the hydrostatic matrix $\boldsymbol {H}$ depends on the local pressure:

$$\nabla\left({\boldsymbol {H}}_k \cdot {\boldsymbol {T}}_v \right)... ...\boldsymbol {q}} \right] + {\boldsymbol {T}}_v \cdot \nabla {\boldsymbol {H}}_k$$ (3.290)

where $\epsilon_v$ is $(R_v/R - 1)$ and $R_v$ is the gas constant for water vapor. The gradient of $T$ is calculated from the spectral representation and that of $q$ from a discrete cubic approximation that is consistent with the interpolation used in the semi-Lagrangian water vapor advection. In general, the elements of $\boldsymbol {H}$ are functions of pressure at adjacent discrete model levels

$$H_{kl} = f_{kl}(p_{l+1/2},p_l,p_{l-1/2})$$ (3.291)

The gradient is then a function of pressure and the pressure gradient

$$\nabla H_{kl} = g_{kl}(p_{_{l+1/2}},\; p_{_{l}},\; p_{_{l-1/2}}, \nabla p_{_{l+1/2}}, \nabla p_{_{l}}, \nabla p_{_{l-1/2}})$$ (3.292)

The pressure gradient is available from (3.263) and the surface pressure gradient calculated from the spectral representation

$$\nabla p_{_{l}} = B_l\nabla p_s = B_l p_s \nabla {\rm ln} p_s$$ (3.293)

3.2.10 Development of semi-implicit system equations

The momentum equation can be written as

$${{\boldsymbol {V}}^{n+1}_{_{A}} - {\boldsymbol {V}}^{n}_{_{D}} \over \Delta t} = -{1 \over 2}\bigg\lbrace \left(1 + \epsilon \right) \left[f{\boldsymbol{\hat{k}}} \times {... ...f{\boldsymbol{\hat{k}}} \times {\boldsymbol {V}} \right]^{n}_{_{D}}\bigg\rbrace$$

$$- {1 \over 2} \bigg\lbrace \left(1 + \epsilon \right) \nabla \left[ R {\boldsymbol {H}}_k^r ... ...RT^r {\rm ln} p_s' \right]^{n+1}_{_{A}}\bigg\rbrace + RHS_{\boldsymbol {V}} ,$$ (3.294)

where $RHS_{\boldsymbol {V}}$ contains known terms at times ( ${n+{1\over 2}}$) and ($n$).

By combining terms, 3.294 can be written in general as

$${\cal U}^{^{n+1}}_{_{A}} {\bf\hat{i}}_{_{A}} + {\cal V}^{^{n+1}}_... ...} + {\cal U}_{_{D}} {\bf\hat{i}}_{_{D}} + {\cal V}_{_{D}} {\bf\hat{j}}_{_{D}} ,$$ (3.295)

where $\bf {\hat{i}}$ and $\bf {\hat{j}}$ denote the spherical unit vectors in the longitudinal and latitudinal directions, respectively, at the points indicated by the subscripts, and $\cal U$ and $\cal V$ denote the appropriate combinations of terms in 3.294. Note that ${\cal U}^{^{n+1}}_{_{A}}$ is distinct from the ${\cal U}_{_{A}}$. Following Bates et al. [13], equations for the individual components are obtained by relating the unit vectors at the departure points ( ${\bf\hat{i}}_{_{D}}$, ${\bf\hat{j}}_{_{D}}$) to those at the arrival points ( ${\bf\hat{i}}_{_{A}}$, ${\bf\hat{j}}_{_{A}}$):

$${\bf\hat{i}}_{_{D}} = \alpha^u_{_{A}}{\bf\hat{i}}_{_{A}} + \beta^u_{_{A}}{\bf\hat{j}}_{_{A}}$$ (3.296)

$${\bf\hat{j}}_{_{D}} = \alpha^v_{_{A}}{\bf\hat{i}}_{_{A}} + \beta^v_{_{A}}{\bf\hat{j}}_{_{A}}  ,$$ (3.297)

in which the vertical components ( ${\bf\hat{k}}$) are ignored. The dependence of $\alpha$'s and $\beta$'s on the latitudes and longitudes of the arrival and departure points is given in the Appendix of Bates et al. [13].

W&O94 followed Bates et al. [13] which ignored rotating the vector to remain parallel to the earth's surface during translation. We include that factor by keeping the length of the vector written in terms of $\left({\boldsymbol{\hat{i}}}_{_{A}},{\boldsymbol{\hat{j}}}_{_{A}}\right)$ the same as the length of the vector written in terms of $\left({\boldsymbol{\hat{i}}}_{_{D}},{\boldsymbol{\hat{j}}}_{_{D}}\right)$. Thus, (10) of W&O94 becomes

$${\cal U}^{^{n+1}}_{_{A}} = {\cal U}_{_{A}} + \gamma\alpha^u_{_{A}}{\cal U}_{_{D}} + \gamma\alpha^v_{_{A}}{\cal V}_{_{D}}$$

$${\cal V}^{^{n+1}}_{_{A}} = {\cal V}_{_{A}} + \gamma\beta^u_{_{A}}{\cal U}_{_{D}} + \gamma\beta^v_{_{A}}{\cal V}_{_{D}}$$ (3.298)

where

$$\gamma = \left[ { {\cal U}_{_{D}}^2 + {\cal V}_{_{D}}^2 \over \le... ...}}\beta^u_{_{A}} + {\cal V}_{_{D}}\beta^v_{_{A}}\right)^2 } \right]^{1 \over 2}$$ (3.299)

After the momentum equation is written in a common set of unit vectors

$${\boldsymbol {V}}^{n+1}_{_{A}} + \left({1 + \epsilon \over 2}\rig... ...{T}} + RT^r {\rm ln} p_s' \right]^{n+1}_{_{A}} = {\cal R}_{\boldsymbol {V}}^*$$ (3.300)

Drop the $(   )^{n+1}_A$ from the notation, define

$$\alpha = \left(1 + \epsilon\right) \Delta t \Omega$$ (3.301)

and transform to vorticity and divergence

$$\zeta + \alpha \sin \varphi \delta + {\alpha \over a} v \cos\varphi = {1 \over a\cos\varphi} \left[ {\partial{\cal R}_v^* \over \partia... ... - {\partial \over \partial\varphi}\left({\cal R}_u^* \cos\varphi\right)\right]$$ (3.302)

$$\delta - \alpha \sin \varphi \zeta + {\alpha \over a} u \cos\varphi + \left( {1 + \epsilon \over 2} \right) \Delta t \nabla^2 \left[ R ... ...ol {H}}_k^r \cdot {\boldsymbol {T}} + RT^r {\rm ln} p_s' \right]^{n+1}_{_{A}}$$

$$= {1 \over a\cos\varphi} \left[ {\partial{\cal R}_u^* \over \partia... ... + {\partial \over \partial\varphi}\left({\cal R}_v^* \cos\varphi\right)\right]$$ (3.303)

Note that

$$u\cos\varphi = {1 \over a} {\partial \over\partial \lambda}\left(\nabla^{-2} \de... ...varphi \over a} {\partial \over\partial \varphi}\left(\nabla^{-2} \zeta \right)$$ (3.304)

$$v\cos\varphi = {1 \over a} {\partial \over\partial \lambda}\left(\nabla^{-2} \ze... ...arphi \over a} {\partial \over\partial \varphi}\left(\nabla^{-2} \delta \right)$$ (3.305)

Then the vorticity and divergence equations become

$$\zeta + \alpha \sin \varphi \delta + {\alpha \over a^2}{\partial \over\partial \lambda}\left(\nabla^{-2} \zeta \right) + {\alpha\cos\varphi \over a^2}{\partial \over\partial \varphi}\left(\nabla^{-2} \delta \right)$$

$$= {1 \over a\cos\varphi} \left[ {\partial {\cal R}_v^* \over \parti... ... \over \partial \varphi}\left({\cal R}_u^* \cos\varphi\right)\right] = {\cal L}$$ (3.306)

$$\delta - \alpha \sin \varphi \zeta + {\alpha \over a^2}{\partial \over\partial \lambda}\left(\nabla^{-2} \delta \right) - {\alpha\cos\varphi \over a^2}{\partial \over\partial \varphi}\lef... ...ol {H}}_k^r \cdot {\boldsymbol {T}} + RT^r {\rm ln} p_s' \right]^{n+1}_{_{A}}$$

$$= {1 \over a\cos\varphi} \left[ {\partial {\cal R}_u^* \over \parti... ... \over \partial \varphi}\left({\cal R}_v^* \cos\varphi\right)\right] = {\cal M}$$ (3.307)

Transform to spectral space as described in the description of the Eulerian spectral transform dynamical core. Note, from (4.5b) and (4.6) on page 177 of Machenhauer [118]

$$\mu P_n^m = D_{n+1}^m P_{n+1}^m + D_{n}^m P_{n-1}^m$$ (3.308)

$$D_{n}^m = \left({n^2-m^2 \over 4n^2-1}\right)^{1 \over 2}$$ (3.309)

and from (4.5a) on page 177 of Machenhauer [118]

$$\left( 1 - \mu^2\right) {\partial \over\partial \mu}P_n^m = -nD_{n+1}^m P_{n+1}^m + \left (n+1 \right) D_{n}^m P_{n-1}^m$$ (3.310)

Then the equations for the spectral coefficients at time $n+1$ at each vertical level are

$$\zeta_n^m \left( 1 - {im\alpha \over n\left(n+1\right)}\right) + ... ...n+1}\right) D_{n+1}^m + \delta_{n-1}^m \alpha \left({n+1\over n}\right) D_{n}^m = {\cal L}_n^m$$ (3.311)

$$\delta_n^m \left( 1 - {im\alpha \over n\left(n+1\right)}\right) -... ... n+1}\right) D_{n+1}^m - \zeta_{n-1}^m \alpha \left({n+1\over n}\right) D_{n}^m$$ (3.312)

$$- \left({1 + \epsilon \over 2}\right) \Delta t {n\left(n+1\right)... ...ymbol {H}}_k^r \cdot {\boldsymbol {T}_n^m} + RT^r {\rm ln} {p_s'}_n^m \right] = {\cal M}_n^m$$

$${{\rm ln} p'_s}^m_n = {\rm PS}^m_n - \left({ 1 + \epsilon \over 2}\right) {\Delta t\over p^r_s} \left(\underline{\Delta p^r} \right)^T \underline{\delta}^m_n$$ (3.313)

$$\underline{T}^m_n = \underline{\rm TS}^m_n - \left({ 1 + \epsilon \over 2}\right) \Delta t {\boldsymbol {D}}^r \underline{\delta}^m_n$$ (3.314)

The underbar denotes a vector over vertical levels. Rewrite the vorticity and divergence equations in terms of vectors over vertical levels.

$$\underline{\delta}_n^m \left( 1 - {im\alpha \over n\left(n+1\righ... ... D_{n+1}^m \underline{\zeta}_{n-1}^m \alpha \left({n+1\over - n}\right) D_{n}^m$$ (3.315)

$$- \left({1 + \epsilon \over 2}\right) \Delta t {n\left(n+1\right)... ...ymbol {H}}^r \underline{T}_n^m + R\underline{T}^r {\rm ln} {p_s'}_n^m \right] = \underline{DS}_n^m$$

$$\underline{\zeta}_n^m \left( 1 - {im\alpha \over n\left(n+1\right... ...D_{n+1}^m + \underline{\delta}_{n-1}^m \alpha \left({n+1\over n}\right) D_{n}^m = \underline{VS}_n^m$$ (3.316)

Define $\underline{h}_n^m$ by

$$g\underline{h}_n^m = R {\boldsymbol {H}}^r T_n^m + R\underline{T}^r {\rm ln} {p_s'}_n^m$$ (3.317)
and

$${\cal A}_n^m = 1 - {im\alpha \over n\left(n+1\right)}$$ (3.318)

$${{\cal B}^+}_n^m = \alpha\left({n\over n+1}\right) D_{n+1}^m$$ (3.319)

$${{\cal B}^-}_n^m = \alpha\left({n+1\over n}\right) D_{n}^m$$ (3.320)

Then the vorticity and divergence equations are

$${\cal A}_n^m \underline{\zeta}_n^m + {{\cal B}^+}_n^m \underline{\delta}_{n+1}^m + {{\cal B}^-}_n^m \underline{\delta}_{n-1}^m = \underline{\hbox{\sffamily\slshape V}\hbox{\sffamily\slshape S}}_n^m$$ (3.321)

$${\cal A}_n^m \underline{\delta}_n^m - {{\cal B}^+}_n^m \underline... ...silon \over 2}\right) \Delta t {n\left(n+1\right) \over a^2} g\underline{h}_n^m = \underline{\hbox{\sffamily\slshape D}\hbox{\sffamily\slshape S}}_n^m$$ (3.322)

Note that these equations are uncoupled in the vertical, i.e. each vertical level involves variables at that level only. The equation for $\underline{h}_n^m$ however couples all levels.

$$g\underline{h}_n^m = -\left({1 + \epsilon \over 2}\right) \Delta ... ...amily\slshape T}\hbox{\sffamily\slshape S}}_n^m + R\underline{T}^r {\rm PS}_n^m$$ (3.323)

Define ${\boldsymbol {C}}^r$ and $\underline{\hbox{\sffamily\slshape H}\hbox{\sffamily\slshape S}}_n^m$ so that

$$g\underline{h}_n^m = -\left({1 + \epsilon \over 2}\right) \Delta ... ...lta}_n^m + \underline{\hbox{\sffamily\slshape H}\hbox{\sffamily\slshape S}}_n^m$$ (3.324)

Let $g{\rm D}_\ell$ denote the eigenvalues of ${\boldsymbol {C}}^r$ with corresponding eigenvectors $\underline{\Phi}_\ell$ and $\boldsymbol{\Phi}$ is the matrix with columns $\underline{\Phi}_\ell$

$$\boldsymbol{\Phi}= \left( \begin{array}{*{3}{c@{\:}}c} \underline{\Phi}_1 & \underline{\Phi}_2 & \dots & \underline{\Phi}_L \ \end{array}\right)$$ (3.325)

and $g{\boldsymbol {D}}$ the diagonal matrix of corresponding eigenvalues

$$g{\boldsymbol {D}} = g \left( \begin{array}{*{3}{c@{\:}}c} {\rm D}_1 & 0 & \cdots ... ...s & \vdots \ 0 & 0 & \cdots & {\rm D}_L \\ \end{array}\right)$$ (3.326)

$${\boldsymbol {C}}^r {\boldsymbol{\Phi}} = {\boldsymbol{\Phi}} g {\boldsymbol {D}}$$ (3.327)

$${\boldsymbol{\Phi}}^{-1}{\boldsymbol {C}}^r {\boldsymbol{\Phi}} = g{\boldsymbol {D}}$$ (3.328)

Then transform

$$\underline{\tilde\zeta}_n^m = \boldsymbol{\Phi}^{-1} \underline{\zeta}_n^m \quad,\quad \underline{\widetilde{\hbox{\sffamily\slshape V}\hbox{\sffamily\slshape S}}}_n^m = \boldsymbol{\Phi}^{-1} \underline{\hbox{\sffamily\slshape V}\hbox{\sffamily\slshape S}}_n^m$$ (3.329)

$$\underline{\tilde\delta}_n^m = \boldsymbol{\Phi}^{-1} \underline{\delta}_n^m \quad,\quad \underline{\widetilde{\hbox{\sffamily\slshape D}\hbox{\sffamily\slshape S}}}_n^m = \boldsymbol{\Phi}^{-1} \underline{\hbox{\sffamily\slshape D}\hbox{\sffamily\slshape S}}_n^m$$ (3.330)

$$\underline{\tilde h}_n^m = \boldsymbol{\Phi}^{-1} \underline{ h}_n^m \quad,\quad \underline{\widetilde{\hbox{\sffamily\slshape H}\hbox{\sffamily\slshape S}}}_n^m = \boldsymbol{\Phi}^{-1} \underline{\hbox{\sffamily\slshape H}\hbox{\sffamily\slshape S}}_n^m$$ (3.331)

$${\cal A}_n^m \underline{\tilde\zeta}_n^m + {{\cal B}^+}_n^m \underline{\tilde\delta}_{n+1}^m + {{\cal B}^-}_n^m \underline{\tilde\delta}_{n-1}^m = \underline{\widetilde{\hbox{\sffamily\slshape V}\hbox{\sffamily\slshape S}}}_n^m$$ (3.332)

$${\cal A}_n^m \underline{\tilde\delta}_n^m - {{\cal B}^+}_n^m \und... ...over 2}\right) \Delta t {n\left(n+1\right) \over a^2} g\underline{\tilde h}_n^m = \underline{\widetilde{\hbox{\sffamily\slshape D}\hbox{\sffamily\slshape S}}}_n^m$$ (3.333)

$$g\underline{\tilde h}_n^m + \left({1 + \epsilon \over 2}\right) \... ...\boldsymbol {C}}^r\boldsymbol{\Phi}\boldsymbol{\Phi}^{-1}\underline{\delta}_n^m = \underline{\widetilde{\hbox{\sffamily\slshape H}\hbox{\sffamily\slshape S}}}_n^m$$ (3.334)

$$\underline{\tilde h}_n^m + \left({1 + \epsilon \over 2}\right) \Delta t {\boldsymbol {D}}\underline{\tilde \delta}_n^m = {1 \over g} \underline{\widetilde{\hbox{\sffamily\slshape H}\hbox{\sffamily\slshape S}}}_n^m$$ (3.335)

Since ${\boldsymbol {D}}$ is diagonal, all equations are now uncoupled in the vertical.

For each vertical mode, i.e. element of $(\underline{\tilde{  }})_n^m$, and for each Fourier wavenumber $m$ we have a system of equations in $n$ to solve. In following we drop the Fourier index $m$ and the modal element index $(  )_\ell$ from the notation.

$${\cal A}_n {\tilde\zeta}_n + {{\cal B}^+}_n {\tilde\delta}_{n+1} + {{\cal B}^-}_n {\tilde\delta}_{n-1} = {\widetilde{\hbox{\sffamily\slshape V}\hbox{\sffamily\slshape S}}}_n$$ (3.336)

$${\cal A}_n {\tilde\delta}_n - {{\cal B}^+}_n {\tilde\zeta}_{n+1} ... ...+ \epsilon \over 2}\right) \Delta t {n\left(n+1\right) \over a^2} g{\tilde h}_n = {\widetilde{\hbox{\sffamily\slshape D}\hbox{\sffamily\slshape S}}}_n$$ (3.337)

$${\tilde h}_n + \left({1 + \epsilon \over 2}\right) \Delta t {\rm D}_\ell{\tilde \delta}_n = {1 \over g} {\widetilde{\hbox{\sffamily\slshape H}\hbox{\sffamily\slshape S}}}_n$$ (3.338)

The modal index $(  )_\ell$ was included in the above equation on ${\rm D}$ only as a reminder, but will also be dropped in the following.

Substitute ${\tilde\zeta}$ and ${\tilde h}$ into the ${\tilde\delta}$ equation.

$$\left[ {\cal A}_n + \left( {1 + \epsilon \over 2}\right)^2 \left(... ...cal B}^-}_{n+1} + {{\cal B}^-}_n {\cal A}_{n-1}^{-1} {{\cal B}^+}_{n-1} \right] {\tilde\delta}_{n} \cr + \left( {{\cal B}^+}_n {\cal A}_{n+1}^{-1... ..._{n-1}^{-1} {{\cal B}^-}_{n-1} \right) {\tilde\delta}_{n-2}$$ (3.339)

$$= \widetilde{\hbox{\sffamily\slshape D}\hbox{\sffamily\slshape S}... ...-1}^{-1} \widetilde{\hbox{\sffamily\slshape V}\hbox{\sffamily\slshape S}}_{n-1}$$

which is just two tri-diagonal systems of equations, one for the even and one for the odd $n$'s, and $m \le n \le N$

At the end of the system, the boundary conditions are

$$n = m, \quad{{\cal B}^-}_n = {{\cal B}^-}_m^m = 0$$ (3.340)

$$n = m+1, \quad {{\cal B}^-}_{n-1} = {{\cal B}^-}_m^m = {{\cal B}^-}_{(m+1)-1}^m = 0$$

the ${\tilde\delta}_{n-2}$ term is not present, and from the underlying truncation

$$\tilde\delta_{N+1}^m = \tilde\delta_{N+2}^m = 0$$ (3.341)

For each $m$ and $\ell$ we have the general systems of equations

$$-A_n {\tilde\delta}_{n+2} + B_n {\tilde\delta}_{n} -C_n -{\tilde\delta}_{n-2} = D_n \;, \left\{ \begin{array}{c} n=m,m+2,..., \left\{ \begin{... ... {\rm or} \ N+2 \\ \end{array}\right. \\ \end{array}\right.$$ (3.342)

$$C_m = C_{m+1} =$$ 0 (3.343)

$$\tilde\delta_{N+1} = \tilde\delta_{N+2} =$$ 0 (3.344)

Assume solutions of the form

$$\tilde\delta_n = E_n \tilde\delta_{n+2} + F_n$$ (3.345)

then

$$E_m = {A_m \over B_m}$$ (3.346)

$$F_M = {D_m \over B_m}$$ (3.347)

$$E_n = {A_n \over B_n - C_nE_{n-2}}  ,    n=m+2,m+4,..., \left\{ \begin{array}{c} N-2 \ {\rm or} \ N-3 \\ \end{array}\right.$$ (3.348)

$$F_n = {D_n + C_nF_{n-2} \over B_n - C_nE_{n-2}}  ,    n=m+2,m+4,...... ...\{ \begin{array}{c} N \ {\rm or} \ N-1 \\ \end{array}\right.$$ (3.349)

$$\tilde\delta_N = F_N    {\rm or}    \tilde\delta_{N-1} = F_{N-1} \;,$$ (3.350)

$$\tilde\delta_n = E_n\tilde\delta_{n+2} + F_n \;,   \left\{ \begin... ...in{array}{c} m+1 \ {\rm or} \ m \ \end{array} \right. \ \end{array} \right.$$ (3.351)

Divergence in physical space is obtained from the vertical mode coefficients by

$$\underline{\delta}_n^m = \boldsymbol{\Phi}\underline{\tilde\delta}_n^m$$ (3.352)

The remaining variables are obtained in physical space by

$$\zeta_n^m \left( 1 - {im\alpha \over n\left(n+1\right)}\right) = {\cal L}_n^m - \delta_{n+1}^m \alpha \left({n\over n+1}\right) D_{n+1}^m - \delta_{n-1}^m \alpha \left({n+1\over n}\right) D_{n}^m$$ (3.353)

$$\underline{T}^m_n = \underline{\rm TS}^m_n - \left({ 1 + \epsilon \over 2}\right) \Delta t {\boldsymbol {D}}^r \underline{\delta}^m_n$$ (3.354)

$${{\rm ln} p'_s}^m_n = {\rm PS}^m_n - \left({ 1 + \epsilon \over 2}\right) {\Delta t\over p^r_s} \left(\underline{\Delta p^r} \right)^T \underline{\delta}^m_n$$ (3.355)

3.2.11 Trajectory Calculation

The trajectory calculation follows Hortal [77] Let ${\boldsymbol {R}}$ denote the position vector of the parcel,

$${d{\boldsymbol {R}} \over dt} = {\boldsymbol {V}}$$ (3.356)

which can be approximated in general by

$${\boldsymbol {R}}_D^n = {\boldsymbol {R}}_A^{n+1} - \Delta t {\boldsymbol {V}}_M^{n+{1 \over 2}}$$ (3.357)

Hortal's method is based on a Taylor's series expansion

$${\boldsymbol {R}}_A^{n+1} = {\boldsymbol {R}}_D^n + \Delta t \lef... ...a t^2 \over 2} \left( {d^2 {\boldsymbol {R}} \over d t^2} \right)_D^{n} + \dots$$ (3.358)

or substituting for ${d {\boldsymbol {R}} / d t}$

$${\boldsymbol {R}}_A^{n+1} = {\boldsymbol {R}}_D^n + \Delta t {\bo... ...Delta t^2 \over 2} \left( {d {\boldsymbol {V}} \over d t} \right)_D^{n} + \dots$$ (3.359)

Approximate

$$\left( {d {\boldsymbol {V}} \over d t} \right)_D^{n} \approx {{\boldsymbol {V}}_A^n - {\boldsymbol {V}}_D^{n-1} \over \Delta t}$$ (3.360)
giving

$${\boldsymbol {V}}_M^{n+{1 \over 2}} = {1 \over 2}\left[\left(2{\boldsymbol {V}}^n-{\boldsymbol {V}}^{n-1}\right)_D + {\boldsymbol {V}}_A^n\right]$$ (3.361)

for the trajectory equation.

3.2.12 Mass and energy fixers and statistics calculations

The semi-Lagrangian dynamical core applies the same mass and energy fixers and statistical calculations as the Eulerian dynamical core. These are described in sections 3.1.19, 3.1.20, and 3.1.21.