Time splitting of the RK steps

Still more CPUs can be economized by time-integrating those terms related to gravity waves using a split time step together with a forward-backward time differencing scheme. The method here is fairly standard. The basic idea is that a split time step is defined such that the small time step respects the stability criteria for gravity waves (hydrostatic mode, acoustic waves when non-hydrostatic dynamics are turned on) while the dynamics time step (for advection for example) can be larger (it respects some criteria related to the advective Courant number). CPUs are economized because only the terms involving gravity or acoustic waves are integrated using the small time step, while the more expensive terms (such as horizontal advection) are integrated using the larger time step. When the non-hydrostatic option is active, ASP uses the so-called HE-VI time integration method (Horizontally Explicit-Vertically Implicit) to solve several equations describing the non-hydrostatic dynamics.

The set of time split dry-mass coordinate prognostic equations including moisture effects can be expressed as

$$\frac{U^{\tau+\Delta\tau}-U^{\tau}}{\Delta\tau} + {PG_u}^{\tau+\Delta\tau/2} = F_u^{t^\ast} + F_u^t$$ (207)

$$\frac{V^{\tau+\Delta\tau}-V^{\tau}}{\Delta\tau} + {PG_v}^{\tau+\Delta\tau/2} = F_v^{t^\ast} + F_v^t$$ (208)

$${\frac{\pi_{d s}^{\tau+\Delta\tau}-\pi_{d s}^\tau}{\Delta\tau}} + \int_0^1 \nabla\cdot {\left(\mu_d{\bf V}\right)}^{\tau+\Delta\tau} \, d\eta = 0$$ (209)

$${\frac{\theta_\rho^{\tau+\Delta\tau}\mu_d^{\tau+\Delta\tau} -\the... ... \left[ {\left(\mu_d{\bf V}\right)}^{\tau+\Delta\tau}\theta_\rho^{\tau} \right] = F_{\theta_\rho}^t$$ (210)

$${\frac{w^{\tau+\Delta\tau}\mu^{\tau+\Delta\tau} -w^{\tau}\mu^{\ta... ...ial{\pi_{\pi d}}^{\tau+\Delta\tau}} \,-\, 1 \right) \,-\, F_w^{\tau+\Delta\tau} = F_w^{t^\ast} + F_w^t$$ (211)

$$\frac{ \Phi^{\tau+\Delta\tau} -\Phi^{\tau}}{\Delta\tau} + \frac{\... ...+\Delta\tau}} \frac{\partial\Phi^\tau}{\partial\eta} - g \, w^{\tau+\Delta\tau} = F_\Phi^{t^\ast}$$ (212)

$$\frac{ r^{\tau+\Delta\tau} {\mu_d}^{\tau+\Delta\tau} - r^{\tau}{\mu_d}^{\tau}}{\Delta\tau} = F_r^{t^\ast}$$ (213)

The pressure gradient terms in Eq.s 209-210 are forward-weighted using the Adams-Bashforth explicit time differencing scheme

$$PG^{\tau+\Delta\tau/2} = \mu_d^\tau \left[ \omega_{PG}\,\chi^\tau + \left(1-\omega_{PG}\right)\,\chi^{\tau-\Delta\tau} \right]$$ (214)

where

$$\chi^\tau = \left(1+\epsilon^\tau\right)\,\nabla\Phi^\tau \,+\, C_{pd}\,\theta_\rho^\tau\,\nabla\Pi^\tau$$ (215)

and $\omega_{PG}=1.5$ (Adams-Bashforth weight).

$r$ is a single variable which represents the sum of all water scalar variables defined in Eq. 43: the prognostic updates for each mixing ratio component is done after the call to the dynamic core.

Eq.s 213-215 correspond to the non-hydrostatic add-on module and they are not evaluated if the non-hydrostatic option is not activated (in this case, $w=0$, $p=\pi$, $\Phi$ is diagnosed from the hydrostatic relationship, and $r$ corresponds to the current RK step value which is constant over the time split). When the non-hydrostatic option is activated, these three equations (together with the diagnostic relationships defined in Eq.s 238-240) are solved simultaneously for $w^{\tau+\Delta\tau}$, $\Phi^{\tau+\Delta\tau}$ and $p^{\tau+\Delta\tau}$ using the HEVI (horizontally-explicit, vertically-implicit) method (i.e. simply using a tri-diagonal matrix approach each split time step). Finally, $w$ has a damping term evaluated at $\tau+\Delta\tau$ which is defined from Eq. 120 as

$$F_w^{\tau+\Delta\tau} = -\kappa_w \,\mu_d^{\tau+\Delta\tau}\,w^{\tau+\Delta\tau}$$ (216)

where $\kappa_w$ only varies in the vertical and is zero except in the uppermost layers where it increases approaching a maximum value (unity) at the upper boundary. This is akin to a Raleigh damping to the basic state value (zero) to avoid reflections off the model top (downward).

The terms on the LHS of Eq.s 209-215 are evaluated for each split time step, $\tau$, and the forcing terms, $F$, on the RHS of these equations are held constant over the time splitting interval. The $F\left(t^\ast\right)$ terms are updated at each RK time step (horizontal advection, vertical wind advection for the $U$, $V$ equations and $\Phi$ equation, horizontal and vertical advection for the $w$ equation). The forcing terms, $F(t)$, are only updated once per large (model) time step (horizontal diffusion, filters, physics lateral boundary tendencies). Note that two physics tendencies are not included in the time-split or RK forcings: they are evaluated after the RK steps are done (large scale precipitation as a super-saturation adjustment), and vertical turbulent diffusion/surface interactions (fully implicit solution in the vertical with dynamics and physics tendencies as source terms). The time-split averaged tendencies over the current RK step are used to advance $\pi_{ds}$, $\theta_\rho$, $u$, $v$, ($w$ and $\Phi$ if the non-hydrostatic option is activated, which are used to diagnose updated values of $p$, $\Pi$ and $\alpha_d$), respectively, for each RK step.

The forcing functions in Eq.s 209-215 are defined as

$$F_w^{t^\ast} = -\nabla\cdot {\left({\bf V} \,\mu_d\,w\right)}^{t^\ast} - \frac... ...tial}{\partial\eta} {\left({\dot\pi} \,w \right)}^{t^\ast} + f_{w,crv}^{t^\ast}$$ (217)

$$F_w^t = \mu_d^t D_{w}^t$$ (218)

$$F_u^{t^\ast} = -\nabla\cdot {\left({\bf V} \,\mu_d \,u\right)}^{t^\ast} - \fra... ...tial}{\partial\eta} {\left({\dot\pi} \,u \right)}^{t^\ast} + f_{u,crv}^{t^\ast}$$ (219)

$$F_u^t = \mu_d^t \left( D_{u}^t \,+\, {\dot{S}}_{u}^t + {\cal{D}}_{h,u}^t \right)$$ (220)

$$F_v^{t^\ast} = -\nabla\cdot {\left({\bf V} \,\mu_d \,v\right)}^{t^\ast} - \fra... ...tial}{\partial\eta} {\left({\dot\pi} \,v \right)}^{t^\ast} + f_{v,crv}^{t^\ast}$$ (221)

$$F_v^t = \mu_d^t \left( D_{v}^t \,+\, {\dot{S}}_{v}^t + {\cal{D}}_{h,v}^t \right)$$ (222)

$$F_{\theta_\rho}^t = \mu_d^t \left(D_{\theta_\rho}^t \,+\, {\dot{S}}_{\theta_\rho}^t + {\cal{D}}_{h,T}^t \right)$$ (223)

$$F_r^{t^\ast} = -\sum_{j=1}^{N_r} \left[ \nabla\cdot {\left({\bf V} \,\mu_d \,r... ... \frac {\partial}{\partial\eta} {\left({\dot\pi}\, r_j\right)}^{t^\ast} \right]$$ (224)

$$F_r^t = \mu_d^t \left[ \sum_{j=1}^{N_r} \left( D_{r_j}^t \,+\, {\dot{S}}_{r_j}^t \right) + {\cal{D}}_{h,r_q}^t \right]$$ (225)

$$F_\Phi^{t^\ast} = \Phi^{t^\ast} {\left(\nabla\cdot {\bf V}\right)}^{t^\ast} \,-\, \nabla\cdot {\left({\bf V}\,\Phi \right)}^{t^\ast}$$ (226)

Note that the scalar mass fluxes, ${\bf V} \,\mu_d \,r_j$ and ${\dot\pi}\, r_j$, have been adjusted to conserve mass of each scalar, $\mu_d \,r_j$, at the end of each RK step. The forcing term in Eq. 214 consists in the horizontal advection term expressed in terms of a divergence and a flux form advection term (as shown in Eq. 228) since it permits the use of the same advection scheme as is used for the other flux-form prognostic variables.

Note that horizontal advection in Eq. 212 is computed at each split step owing to the strong dependence on the horizontal divergence. When running high resolution runs, the number of time splits can increase significantly depending on the time step size. Thus there is an option to save CPUs when using open lateral boundary conditions (such as in LAM mode): Eq. 212 can be expressed as

$${\frac{\theta_\rho^{\tau+\Delta\tau}\mu_d^{\tau+\Delta\tau} -\the... ...{\bf V}\right)}^{\tau+\Delta\tau} = F_{\theta_\rho}^t + F_{\theta_\rho}^{t\ast}$$ (227)

where the forcing function which is updated at each RK step is expressed as

$$F_{\theta_\rho}^{t\ast} = -\nabla\cdot \left[ {\left(\mu_d{\bf V}... ...\right] + \theta_\rho^{t\ast}\, \nabla\cdot {\left(\mu_d{\bf V}\right)}^{t\ast}$$ (228)

thus the flux form advection (first term on the RHS of Eq. 230) is only evaluated once per RK step (as is the case for all of the other variables). To maintain numerical stability, only the divergence component of the horizontal advection must be updated on the split step, thus it is removed at the RK step and added back at each split step. This form can save CPUs but is essentially only useful for very high (sub-kilometric) spatial resolutions since even for kilometric resolutions, the number of splits is generally at or close to the minimum (default 5 per large time step, $\Delta t$). This is particularly useful for saving CPUs if the option to compute the more computationally intensive 5th order WENO horizontal advection for $\theta_\rho$ at the last RK step only has been activated. Note that the use of Eq. 229 in place of Eq. 212 no longer conserves $\theta_\rho\,\mu_d$, but in LAM mode this is already the case owing to the open lateral boundary source/sink terms. Numerical experiments indicate that in practice, use of Eq. 229 has a negligible impact on the results, even for high resolution runs.

The following diagnostics at time $\tau+\Delta\tau$ used for the solution of the above equations can be defined:

$$\mu_d^{\tau+\Delta\tau} = {\frac{\partial A_\eta}{\partial\eta}}+ {\frac{\partial B_\eta}{\partial\eta}} {\pi_{d s}}^{\tau+\Delta\tau}$$ (229)

$$\pi_{d}^{\tau+\Delta\tau} = A_\eta + B_\eta {\pi_{d s}}^{\tau+\Delta\tau}$$ (230)

$$\Pi_{\pi d}^{\tau+\Delta\tau} = {\left(\frac{\pi_d^{\tau+\Delta\tau}}{p_0}\right)}^{R_d/C_{pd}}$$ (231)

$$\nabla\cdot {\left(\mu_d{\bf V}\right)}^{\tau+\Delta\tau} = \frac{\partial U^{\tau+\Delta\tau}}{\partial x} + \frac{\partial V^{\tau+\Delta\tau}}{\partial y}$$ (232)

$${\dot\pi}^{\tau+\Delta\tau} = B_\eta \int_0^1 {\left( \nabla \cdot {\bf V}\mu_d \right)}^{\ta... ...int_0^\eta {\left( \nabla \cdot {\bf V}\mu_d \right)}^{\tau+\Delta\tau} \,d\eta$$ (233)

$$u^{\tau+\Delta\tau} = {\frac {U^{\tau+\Delta\tau}}{\mu_d^{\tau+\Delta\tau}}}$$ (234)

$$v^{\tau+\Delta\tau} = {\frac {V^{\tau+\Delta\tau}}{\mu_d^{\tau+\Delta\tau}}}$$ (235)

$$p^{\tau+\Delta\tau} = p^{t^\ast} \left[ 1 \,+\, \frac{C_{vd}}{C_{pd}} \left( \frac{\t... ...^{t^\ast}} - \frac{1+r^{t^\ast}}{1+r^{\tau+\Delta\tau}} \,+\, 1 \right) \right]$$ (236)

$$\alpha_d^{\tau+\Delta\tau} = \frac{-1}{{\vartheta_{\pi\,d}}^{\tau+\Delta\tau}} \frac{\partial\Phi^{\tau+\Delta\tau}} {\partial{\Pi_{\pi\,d}}^{\tau+\Delta\tau}}$$ (237)

$$\varphi_q^{\tau+\Delta\tau} = \frac{1}{1 + r^{\tau+\Delta\tau}}$$ (238)

$$\Pi^{\tau+\Delta\tau} = \Pi^{t^\ast} \,+\, \varsigma \frac{\Pi^{t^\ast}}{p^{t^\ast}} \left( p^{\tau+\Delta\tau} - p^{t^\ast} \right)$$ (239)

Eq. 238 is based on the linearization of the ideal gas law (Eq. 37) about the latest Runge-Kutta value for the split time-step integration, the hydrostatic relationship in Eq. 239, and Eq. 240 is the dry to moist density ratio defined in Eq. 43. The diagnostic non-hydrostatic Exner from Eq. 241 is based on a linearization about the latest Runge-Kutta value for the split time-step integration as a function of the non-hydrostatic pressure, $p$. Eq. 37, about the latest Runge-Kutta value for the split time-step integration. $t$ indicates the value at the start of the large current time step, $\Delta t$, and $t^{\ast}$ represents a value at the the current RK step.

The mixing ratios of the different water species, $r$, and other tracers (such as the TKE, $e$), are integrated using the RK tendencies. Note that the use of Eq. 215 implies that the total water mass scalar, $r\,\mu_d$ value is advanced within the same RK time step before the dynamics call, so it uses the same mass fluxes (from the previous RK step) as for the horizontal advection of $\Phi$ and both horizontal and vertical mass fluxes as $u$, $v$ and $w$. Thus, Eq. 215 not only takes into account water mass changes over the time split, but also changes in hydrostatic pressure to compute $\varphi_\rho^{\tau+\Delta\tau}$ so that it is more than just a simple linear interpolation of $r$ over the time split. For large scale or even mesoscale motions, the pertberbation of water vapor or total water scalars can likely be neglected during the time splitting (it is sufficient to update the moisture metrics at each RK step), but for very high resolution runs (such as large-eddy simulations: sub kilometric scales), the integration of $r_v$ or $r$ within the time splitting for computing $p$ becomes more critical (which is the justification for using a water vapor dependent potential temperature, $\theta_m$, in WRF version 4 verses previous versions using $\theta$ as the prognostic variable for temperature). Note that numerical tests with ASP have shown that computing the water vapor scalar RK tendencies using fluxes from the start of the current RK step (verses after the dynamics call, thus using fluxes at the end of the current RK step), has virtually no impact on the results for either hydrostatic or non-hydrostatic simulations, even those at very high resolution. More importantly (to justify the use of computing scalar updates during the time split), any impact has been found to be quite small compared to the impact of using Eq. 215, which is essential, for very high resolution runs.

When the model is run in hydrostatic mode, the overall governing equations are similar except that $\epsilon=0$ in Eq.s 30-31, and Eq.s 32 and Eq. 35 are no longer solved (i.e. $w$ and $\phi_d$ are no longer prognostic). In this case, $w$ can be set to 0 or diagnosed from Eq. 35, but this is of little importance since it is not used in any of the dynamics computations in this case. The pressure is simply diagnosed as $p=\pi$ consistent with $\epsilon=0$, where $\pi$ is the moist hydrostatic pressure (computed from the moist specific density, Eq. 29).The dry inverse density (used to compute the geopotentials) is then simply computed from the ideal gas law (Eq. 37) as

$$\alpha_d = \frac{\theta_\rho\Pi\,R_d}{p\,\varphi_q}$$ (240)

where it follows that for a hydrostatic atmosphere, the pressure is defined as

$$\frac{\partial p}{\partial\eta} = \frac{\mu_d}{\varphi_q} = \mu$$ (241)

Eq. 243 is also used to initialize $p$ when running in non-hydrostatic mode (and correspondingly $w=0$ initially). Finally, map factors are taken into account in the usual manner, along with full curvature terms in the horizontal momentum equations and this is detailed in a subsequent section.

Eq.s 213-239 can be manipulated in finite difference form to yield a set of three equations with three unknowns ( ${\Phi}^{\tau+\Delta\tau}$, $w^{\tau+\Delta\tau}$ and $p^{\tau+\Delta\tau}$):

$$\begin{align*}\begin{split} w_k^{\tau+\Delta\tau} =& \Biggr\lbrace w_k^{\tau} \l... ... \right] \Biggr\rbrace \frac{1}{\left(1+\kappa_w\right)} \end{split}\end{align*}$$ (242)

$$\begin{align*}\begin{split} {\alpha_{d,k}}^{\tau+\Delta\tau}=& \Biggr\lbrace {\h... ...tau+\Delta\tau}- \Pi_{\pi,d,k}^{\tau+\Delta\tau}\right)} \end{split}\end{align*}$$ (243)

$$\begin{align*}\begin{split} p_k^{\tau+\Delta\tau} =& p_k(t^\ast) \Bigg\lbrace 1 ... ...i_{q,k}\left(t^\ast\right)} \,+\, 1 \right] \Bigg\rbrace \end{split}\end{align*}$$ (244)

where Eq. 245 is obtained by substitution of the geopotential tendency Eq. 214 into the hydrostatic equation Eq. 239, and the diagnostic pressure equation, Eq. 246 is obtained by substituting the density ratio Eq. 240 into the linearized gas law Eq. 238.

The forcing functions Eq.s 219 -228 are expanded in terms of the vertical finite differencing as

$$\begin{align*}\begin{split} F_{w,k} \left(t^\ast\right) =& -\nabla\cdot \left( {... ...right]} {\left(\Delta\eta_{k}+\Delta\eta_{k-1}\right)/2} \end{split}\end{align*}$$ (245)

$$\begin{align*}\begin{split} F_{w,k} (t) = & {\hat\mu}_{d,k}\,\left( D_{w\,k} + f_{w,crv,k}\right) \end{split}\end{align*}$$ (246)

$$\begin{align*}\begin{split} F_{u,k} (t^\ast) = & -\nabla\cdot \left( {\bf {V}}_k... ...}^{t\ast} {\hat{u}}_{k}^{t\ast} \right) } {\Delta\eta_k} \end{split}\end{align*}$$ (247)

$$\begin{align*}\begin{split} F_{u,k} (t) =& \mu_{d,k}\left(D_{u,k} + {\dot{S}}_{u,k} + {\cal{D}}_{h,u,k} \right) \end{split}\end{align*}$$ (248)

$$\begin{align*}\begin{split} F_{v,k} (t^\ast) = & -\nabla\cdot \left( {\bf {V}}_k... ...}^{t\ast} {\hat{v}}_{k}^{t\ast} \right) } {\Delta\eta_k} \end{split}\end{align*}$$ (249)

$$\begin{align*}\begin{split} F_{v,k} (t) =&\mu_{d,k}\left(\,D_{v,k} + {\dot{S}}_{v,k} + {\cal{D}}_{h,v,k} \right) \end{split}\end{align*}$$ (250)

$$\begin{align*}\begin{split} F_{\theta_\rho,k} (t) = & \mu_{d,k}\left(D_{\theta_\rho\,k} + {\dot{S}}_{\theta,k} + {\cal{D}}_{h,T,k} \right) \end{split}\end{align*}$$ (251)

$$\begin{align*}\begin{split} F_{r,k}^{t^\ast} = & -\sum_{j=1}^{N_r} \left[ \nabla... ...} {\hat{r}}_{k}^{t\ast} \right) } {\Delta\eta_k} \right] \end{split}\end{align*}$$ (252)

$$\begin{align*}\begin{split} F_r^t = & \mu_{d,k}^t \left[ \sum_{j=1}^{N_r} \left(... ...{S}}_{r_j,k}^t + \right) + {\cal{D}}_{h,r_q,k}^t \right] \end{split}\end{align*}$$ (253)

$$\begin{align*}\begin{split} F_{\alpha,k} \left(t^\ast\right) =& F_{\Phi,k}(t^\as... ...t{\bf V}}_{k+1}}^{t\ast} {\hat\Phi}_{k+1}^{t\ast}\right) \end{split}\end{align*}$$ (254)

while the forcing functions, $F_{w}(t)$ and $F_{\alpha}(t)$, potentially include diffusion terms which have been evaluated on the large dynamics time step, $\Delta t$, before the RK steps.

Eq.s 244-246 are solved simultaneously using a tri-diagonal approach. This is called the HEVI (horizontally explicit and vertically implicit) solution method for the non-hydrostatic dynamics. These three equations (Eq.s 244-246) can then be rewritten in a simplified manner in terms of three unknowns as

$${w_k}^{\tau+\Delta\tau} = {\cal F}_{w,k} + {\cal G}_{w,k} \left({p_k}^{\tau+\Delta\tau}-{p_{k-1}}^{\tau+\Delta\tau}\right)$$ (255)

$$\alpha_{d,k}^{\tau+\Delta\tau} = {\cal F}_{\alpha_k} + {\cal G}_{\alpha} \left({w_k}^{\tau+\Delta\tau}-{w_{k+1}}^{\tau+\Delta\tau}\right)$$ (256)

$${p_k}^{\tau+\Delta\tau} = {\cal F}_{p,k} + {\cal G}_{p,k} \,\alpha_{d,k}^{\tau+\Delta\tau}$$ (257)

Eq.s 257-259 are solved from $k=1,N$. The coefficients of these three linear equations are expressed as

$$\begin{align*}\begin{split} {\cal F}_{w,k} =& \Biggr\lbrace {w_k}^{\tau} \left( ... ...ight] - g\,\Delta\tau \Biggr\rbrace \frac{1}{1+\kappa_w} \end{split}\end{align*}$$ (258)

$$\begin{align*}\begin{split} {\cal G}_{w,k} =& \frac{g\,\Delta\tau\,{\hat{\varphi... ...,k-1}}^{\tau+\Delta\tau}\right) \left(1+\kappa_w\right)} \end{split}\end{align*}$$ (259)

$$\begin{align*}\begin{split} {\cal F}_{\alpha_k} =& \Biggr\lbrace {\hat\Phi}_{k}^... ...tau+\Delta\tau}- \Pi_{\pi,d,k}^{\tau+\Delta\tau}\right)} \end{split}\end{align*}$$ (260)

$$\begin{align*}\begin{split} {\cal G}_{\alpha} =& \frac{g\,\Delta\tau}{\vartheta_... ...tau+\Delta\tau}- \Pi_{\pi,d,k}^{\tau+\Delta\tau}\right)} \end{split}\end{align*}$$ (261)

$$\begin{align*}\begin{split} {\cal F}_{p,k} =& p_k(t^\ast) \Bigg\lbrace 1 \,+\, \... ...i_{q,k}\left(t^\ast\right)} \,+\, 1 \right] \Bigg\rbrace \end{split}\end{align*}$$ (262)

$$\begin{align*}\begin{split} {\cal G}_{p,k} =& - \frac{p_k(t^\ast)}{\gamma\,\alpha_{d,k}(t^\ast)} \end{split}\end{align*}$$ (263)

To simultaneously solve Eq.s 257-259, we first derive an equation for $w$. First, $p$ is eliminated from the $w$ equation by substitution of Eq. 259 into Eq. 257 yielding

$${w_k}^{\tau+\Delta\tau} = {\cal F}_{w,k} + {\cal G}_{w,k} \bigg\lbrack {\cal F}_{p,k} - {\c... ...tau+\Delta\tau} - {\cal G}_{p,k-1}\alpha_{d,k-1}^{\tau+\Delta\tau} \bigg\rbrack$$ (264)

Next, $\alpha_d$ is eliminated by substituting Eq. 258 into Eq. 266 to have

$$\begin{align*}\begin{split} {w_k}^{\tau+\Delta\tau} =& {\cal F}_{w,k} + {\cal G}... ...u} - w_{k}^{\tau+\Delta\tau} \right)\right] \bigg\rbrace \end{split}\end{align*}$$ (265)

This equation can then be solved as a tri-diagonal matrix so that

$${w_k}^{\tau+\Delta\tau} = d_k + a_k {w_{k+1}}^{\tau+\Delta\tau} +... ...u+\Delta\tau} + c_k {w_{k-1}}^{\tau+\Delta\tau} \hskip.5in \left( k=2,N \right)$$ (266)

where the matrix coefficients are defined as

$$a_k = - {\cal G}_{w,k} {\cal G}_{p,k} {\cal G}_{\alpha_k}$$ (267)

$$b_k = {\cal G}_{w,k} \left( {\cal G}_{p,k-1}{\cal G}_{\alpha_{k-1}} + {\cal G}_{p,k} {\cal G}_{\alpha_k}\right)$$ (268)

$$c_k = - {\cal G}_{w,k} {\cal G}_{p,k-1} {\cal G}_{\alpha_{k-1}}$$ (269)

$$d_k = {\cal F}_{w,k} + {\cal G}_{w,k} \left( {\cal F}_{p,k} - {\cal F}_... ...al G}_{p,k}{\cal F}_{\alpha_k} -{\cal G}_{p,k-1}{\cal F}_{\alpha_{k-1}} \right)$$ (270)

Upper Boundary Conditions

$w_k$ and the corresponding $\epsilon_k$ (the vertical pressure gradient) are defined at levels. At the upper boundary ($k=1$), it is assumed that ${\hat{p}}_1=\pi_1$ where ${\hat{p}}_1$ is the non-hydrostatic pressure at the top of the model domain (which is not explicitly used by the model), thus, $dw/dt$ is not necessarily zero in ASP. Assuming that the pressure at the top of the model domain is equal to the hydrostatic pressure using a one-sided finite difference, one can write

$$\epsilon_{1} = {\varphi}_1 \left(\frac{p_1 - \pi_{1}}{{\overline\pi}_{d,1} - \pi_{d,1}}\right) - 1$$ (271)

so Eq. 257 can therefore be expressed as

$${w_1}^{\tau+\Delta\tau} = {\cal F}_{w,1} + {\cal G}_{w,1} \,{p_1}^{\tau+\Delta\tau}$$ (272)

recalling that $\pi_{1}$ is constant in time. Equations 258-259 are unchanged. The ${\cal F}$ coefficients for $w$ and $\alpha_d$ of Eq. 274 are modified at $k=1$ and become

$${\cal F}_{w,1} = \Biggr\lbrace {w_1}^{\tau} \left( \frac{ {{\hat\mu}_{d,1}}^{\tau}... ...ime \tau+\Delta t}-\pi_1\right)} + 1 \right] \Biggr\rbrace \frac{1}{1+\kappa_w}$$ (273)

$${\cal F}_{\alpha_1} = % \Biggr\lbrace {\hat\Phi}_{1}^{\tau} - {\hat\Phi}_{2}^{\tau} + \... ...\left( \Pi_{\pi,d,2}^{\tau+\Delta\tau}- \Pi_{\pi,d,1}^{\tau+\Delta\tau}\right)}$$ (274)

where the RK $w$ forcing term (Eq. 247) is

$$F_{w,1} \left(t^\ast\right) = -\nabla\cdot \left( {{\hat{\bf V}}_... ...st} {w_1}^{t\ast} \right) - \frac{ {\dot{\pi}}_2\,w_2^{t\ast}} {\Delta\eta_{1}}$$ (275)

The upper boundary conditions for the tri-diagonal matrix are

$$a_1 = -{\cal G}_{w,1} {\cal G}_{p,1} {\cal G}_{\alpha_1}$$ (276)

$$b_1 = {\cal G}_{w,1} {\cal G}_{p,1} {\cal G}_{\alpha_1}$$ (277)

$$c_1 =$$ 0 (278)

$$d_1 = {\cal F}_{w,1} + {\cal G}_{w,1} \left( {\cal F}_{p,1} +{\cal G}_{p,1}{\cal F}_{\alpha_1} \right)$$ (279)

Lower Boundary Conditions

The lower boundary condition for Eq. 268 is simply $w_{N+1}=0$ ($dw/dt=0$ right at the surface), which implies that $\epsilon=1$ at the bottom of the model domain. The lower boundary condition for the tri-diagonal matrix consists in simply

$$a_N = 0$$ (280)

and $b_N$, $c_N$, and $d_N$ are as defined previously except for $k=N$.

Once Eq. 268 has been solved for ${w_k}^{\tau+\Delta\tau}$, then it can be substituted into Eq. 258 to obtain ${\alpha_{k}}^{\tau+\Delta\tau}$, which can then be substituted into Eq. 259 to obtain ${p_k}^{\tau+\Delta\tau}$. Finally, the geopotential on levels can then be updated from the hydrostatic equation, Eq. 239.

Time split steps

The split time step obeys the constraint

$$\Delta \tau < \frac{d }{\sqrt{2} \, c_{max}}$$ (281)

where $c_{max}$ is the speed of sound. But in fact, this values has no impact on the large dynamics time step but rather controls the number of time splits per RK step (again, requiring the evaluation of only those terms related to gravity or sound wave propagation). The WRF model currently recommends:

$$\Delta \tau = \frac{C_\tau \, d}{2\, c_{max}}$$ (282)

in order to be compatible with non-hydrostatic scales, thus ASP adopts this criteria. In addition, WRF recommends an additional safety margin using the leading coefficient $C_\tau$ which has a value less than 1: the current ASP default in 0.8. This relatively simple splitting method is run using a model (large) time step of

$$\Delta t = \frac{C_{\Delta t}\,d}{\sqrt{3}\,V_{max}}$$ (283)

where $V_{max}$ is a user-defined maximum wind speed expected in the simulation. The leading coefficient $C_{\Delta t}$ here currently uses a value of 0.75 (based on considerations outlined in the WRFv3 and v4 scientific documentation). The time splitting method used herein represents a trade-off between simplicity (in terms of algorithm, code and a minimal number of computations) verses an increase in the large model time step (using a more complex and computationally expensive time splitting method). As an example, the rule of thumb for MM5 users is to use a time step which is 3 $d$ (km). For example, in ASP assuming $d=35$ km and $V_{max}=100$ m s$^{-1}$, the model time step is defined as $\Delta t=4.3\,d$ (km). The WRF model Skamarock et al. (2005) defines, for 6th order horizontal advection, a large model time step of 4.67 $d$ (km). If the 5th order advection scheme is used in WRF, then the large model time step can be roughly 6 $d$ (km). This represents an improvement over MM5, although more computations are required: indeed, owing to a rather elaborate linearization procedure used by WRF there are numerous terms to evaluate. The relatively simple and efficient splitting method outlined herein is used as the time step values fall between those used by MM5 and WRF.