The nonhydrostatic governing equations of the atmosphere on a sphere read
\begin{eqnarray} &&\dfrac{du}{dt}=-\dfrac{c_p\theta}{r\cos\varphi}\dfrac{\partial\Pi}{\partial\lambda}+f_rv-f_\varphi w+\left(\dfrac{uv\tan\varphi}{r} -\dfrac{uw}{r}\right),\\ &&\dfrac{dv}{dt}=-\dfrac{c_p\theta}{r}\dfrac{\partial\Pi}{\partial\varphi}-f_ru+f_\lambda w-\left(\dfrac{u^2\tan\varphi}{r}+\dfrac{vw}{r}\right),\qquad\quad\\ &&\dfrac{dw}{dt}=-c_p\theta\dfrac{\partial\Pi}{\partial r}-g+f_\varphi u-f_\lambda v+\left(\dfrac{u^2\tan\varphi}{r}+\dfrac{vw}{r}\right),\\ &&(\gamma-1)\dfrac{d\Pi}{dt}=-\Pi D_3+\dfrac{F_\theta^*}{\theta} ,\\ &&\dfrac{d\theta}{dt}=\dfrac{F_\theta^*}{\Pi} , \end{eqnarray}
on a spherical coordinate system, where \(\Pi\) is the Exner function, θ is the potential temperature, cp=1004.64 J kg-1 K-1 represents the specific heat at constant pressure, u, v are horizontal winds and w is the vertical speed, t means time, g=9.80616 m s-2 is the gravitational acceleration, r is the radius vector of the spherical coordinate, D3 denotes the 3D divergence, f is Coriolis parameter and
\begin{eqnarray*} \gamma&=&\dfrac{c_p}{R_{\rm d}} ,\\ F_\theta^*&=&\dfrac{Q_{\rm T}+F_{\rm T}}{c_p} . \end{eqnarray*}
R d represents the ideal gas constant for dry air, Q T shows the source or sink term of heat and F T is the turbulent diffusion, both of which are 0 in this dry core.
For uniform code design on the Yin and Yang components, 3D Coriolis force, different from the original GRAPES (Chen et al., 2008), is described in the momentum equations. It also serves to improve the accuracy and flexibility of the dynamic core. After discretization with the SISL method, the dynamical equations for the prognostic variables \(u,v,w,\Pi'\) and θ' at time level n+1 are written as
\begin{eqnarray} \label{eq3} u_{n+1}&=&\alpha_\varepsilon(L_u)_{n+1}\Delta t+A_u ,\\ \label{eq4} v_{n+1}&=&\alpha_\varepsilon(L_v)_{n+1}\Delta t+A_v ,\\ \label{eq5} w_{n+1}&=&\alpha_\varepsilon(L_{\hat{w}})_{n+1}\Delta t+A_{\hat{w}} ,\\ \label{eq6} (\Pi')_{n+1}&=&\alpha_\varepsilon(L_\Pi)_{n+1}\Delta t+A_\Pi ,\\ \label{eq7} (\theta')_{n+1}&=&\alpha_\varepsilon(L_\theta)_{n+1}\Delta t+A_\theta , \end{eqnarray}
in a height-based terrain-following coordinate \(\hatz\), where \(\Pi'\)(θ') is a perturbation of the Exner function (potential temperature) from its reference state \(\tilde\Pi\)(\(\tilde\theta\)). We note that the curvature terms in the momentum equations disappear when a semi-Lagrangian algorithm is used for transporting computation of the 3D vector. The reference state is a height-dependent function of each variable, which is the horizontal average of the given initial fields and is different from that in the original model. In the equations, ∆ t is the time step and αε the contribution adjustment factor. Information on the departure point at time level n and the nonlinear terms are included in Ax(\(x=u,v,w,\Pi',\theta'\)), which remains the same as in the original GRAPES model (Xue and Chen, 2008). Linear terms Lx(\(x=u,v,w,\Pi',\theta'\)) at time level n+1 are
\begin{eqnarray} \label{eq8} L_u&=&-c_p\tilde{\theta}\left(\dfrac{\partial\Pi'}{\partial x}+Z_{\rm sx}\dfrac{\partial\Pi'}{\partial\hat{z}}\right)+f_rv-f_\varphi w ,\\ \label{eq9} L_v&=&-c_p\tilde{\theta}\left(\dfrac{\partial\Pi'}{\partial y}+Z_{\rm sy}\dfrac{\partial\Pi'}{\partial\hat{z}}\right)-f_ru+f_\lambda w ,\\ \label{eq10} L_{\hat{w}}&=&-c_p\tilde{\theta}Z_{\rm st}\dfrac{\partial\Pi'}{\partial\hat{z}}-s-c_p\theta'Z_{\rm st}\dfrac{\partial\tilde{\Pi}}{\partial\hat{z}} +f_\varphi u-f_\lambda v ,\quad\\ \label{eq11} L_\Pi&=&-\dfrac{\partial\tilde{\Pi}}{\partial\hat{z}}\hat{w}-\dfrac{\tilde{\Pi}}{\gamma-1}D_3|_{\hat{z}}-\bigg[\dfrac{z_{\rm T}-\hat{z}}{z_{\rm T}-z_{\rm s}} \dfrac{\partial\tilde{\Pi}}{\partial\hat{z}}-\nonumber\\ &&\dfrac{\tilde{\Pi}}{(\gamma-1)(z_{\rm T}-z_{\rm s})}\bigg](u\phi_{\rm sx}+v\phi_{\rm sy}) ,\\ \label{eq12} L_\theta&=&-w\dfrac{\partial\tilde{\Pi}}{\partial z} , \end{eqnarray}
where
$$ s=\dfrac{z_{\rm T}}{z_{\rm T}-z_{\rm s}}c_p\tilde{\theta}\dfrac{\partial\tilde{\Pi}}{\partial\hat{z}}+g , $$
is a residual fraction of the reference atmosphere from the hydrostatic balance state, and
\begin{eqnarray*} Z_{\rm sx}&=&-\dfrac{z_{\rm T}-\hat{z}}{z_{\rm T}-z_{\rm s}}\phi_{\rm sx} ,\\ Z_{\rm sy}&=&-\dfrac{z_{\rm T}-\hat{z}}{z_{\rm T}-z_{\rm s}}\phi_{\rm sy} ,\\ Z_{\rm st}&=&-\dfrac{z_{\rm T}}{z_{\rm T}-z_{\rm s}} , \end{eqnarray*}
φ sx and φ sy are the topographic slope, z T is the model top, and z s is surface height. The terms containing fφ and fΛ are newly introduced into the equations in comparison to the GRAPES model. Considering Eqs. (9)-(11), the linear equation set of Eqs. (4)-(6) can be solved accordingly:
\begin{eqnarray} \label{eq13} u_{n+1}&=&\xi_{u1}\dfrac{\partial\Pi'}{\partial x}+\xi_{u2}\dfrac{\partial\Pi'}{\partial y}+\xi_{u3}\dfrac{\partial\Pi'}{\partial\hat{z}}+\xi_{u0} ,\\ \label{eq14} v_{n+1}&=&\xi_{v1}\dfrac{\partial\Pi'}{\partial x}+\xi_{v2}\dfrac{\partial\Pi'}{\partial y}+\xi_{v3}\dfrac{\partial\Pi'}{\partial\hat{z}}+\xi_{v0} ,\\ \label{eq15} \hat{w}_{n+1}&=&\xi_{\hat{w}1}\dfrac{\partial\Pi'}{\partial x}+\xi_{\hat{w}2}\dfrac{\partial\Pi'}{\partial y}+\xi_{\hat{w}3} \dfrac{\partial\Pi'}{\partial\hat{z}}+\xi_{\hat{w}0} , \end{eqnarray}
where
\begin{eqnarray} \label{eq16} \xi_{u1}&=&-\dfrac{C_{11}\delta c_p\tilde{\theta}}{C_0} ,\\ \label{eq17} \xi_{u2}&=&-\dfrac{C_{12}\delta c_p\tilde{\theta}}{C_0} ,\\ \label{eq18} \xi_{u3}&=&-\dfrac{(C_{11}Z_{\rm sx}+C_{12}Z_{\rm sy}+C_{13}Z_{\rm st})\delta c_p\tilde{\theta}}{C_0} ,\\ \label{eq19} \xi_{u0}&=&\dfrac{C_{11}A_u+C_{12}A_v+C_{13}A_w-C_{13}\delta\left(s+c_pZ_{\rm st}\frac{\partial\tilde{\Pi}}{\partial\hat{z}}A_\theta\right)}{C_0} ,\nonumber\\\\ \label{eq20} \xi_{v1}&=&-\dfrac{C_{21}\delta c_p\tilde{\theta}}{C_0} , \end{eqnarray} \begin{eqnarray} \label{eq21} \xi_{v2}&=&-\dfrac{C_{22}\delta c_p\tilde{\theta}}{C_0} ,\\ \label{eq22} \xi_{v3}&=&-\dfrac{(C_{21}Z_{\rm sx}+C_{22}Z_{\rm sy}+C_{23}Z_{\rm st})\delta c_p\tilde{\theta}}{C_0} ,\\ \label{eq23} \xi_{v0}&=&\dfrac{C_{21}A_u+C_{22}A_v+C_{23}A_w-C_{23}\delta\left(s+c_pZ_{\rm st}\frac{\partial\tilde{\Pi}}{\partial\hat{z}}A_\theta\right)}{C_0} ,\nonumber\\ \end{eqnarray}
and
\begin{eqnarray} \label{eq24} \xi_{\hat{w}1}&=&-\dfrac{(C_{31}-\alpha_{w1}C_{11}-\alpha_{w2}C_{21})\delta c_p\tilde{\theta}}{\alpha_{w3}C_0} ,\\ \label{eq25} \xi_{\hat{w}2}&=&-\dfrac{(C_{32}-\alpha_{w1}C_{12}-\alpha_{w2}C_{22})\delta c_p\tilde{\theta}}{\alpha_{w3}C_0} ,\\ \label{eq26} \xi_{\hat{w}3}&=&-\bigg[\dfrac{(C_{31}-\alpha_{w1}C_{11}-\alpha_{w2}C_{21})Z_{\rm sx}}{\alpha _{w3} C_0}+\nonumber\\ &&\dfrac{(C_{32}-\alpha_{w1}C_{12}-\alpha_{w2}C_{22})Z_{\rm sy}+}{\alpha _{w3} C_0}\nonumber\\ &&\dfrac{(C_{33}-\alpha_{w1}C_{13}-\alpha_{w2}C_{23})Z_{\rm st}}{\alpha _{w3} C_0}\bigg]\delta c_p \tilde{\theta} ,\\ \label{eq27} \xi_{\hat{w}0}&=&\dfrac{(C_{31}-\alpha_{w1}C_{11}-\alpha_{w2}C_{21})A_u}{\alpha_{w3}C_0}+\nonumber\\ &&\dfrac{(C_{32}-\alpha_{w1}C_{12}-\alpha_{w2}C_{22})A_v}{\alpha_{w3}C_0}+\nonumber\\ &&\dfrac{(C_{33}-\alpha_{w1}C_{13}-\alpha_{w2}C_{23})A_w}{\alpha_{w3}C_0}-\nonumber\\ &&\dfrac{(C_{33}-\alpha_{w1}C_{13}-\alpha_{w2}C_{23})\delta\left(s+c_pZ_{\rm st}\frac{\partial\tilde{\Pi}}{\partial \hat{z}}A_\theta\right)} {\alpha_{w3}C_0} .\quad\ \ \end{eqnarray}
δ=αε∆ t is defined, and matrix C is given as
$$ \bm{C}=\!\!\left(\!\! \begin{array}{c@{\ \ }c@{\ \ }c} \dfrac{\eta+(\delta f_\lambda)^2}{C_0} & \dfrac{\eta\delta f_r+\delta^2f_\lambda f_\varphi}{C_0} & -\dfrac{\delta f_\varphi-\delta^2f_\lambda f_r}{C_0}\\[3mm] \dfrac{\delta^2f_\lambda f_\varphi-\eta\delta f_\varphi}{C_0} & \dfrac{\eta+\eta(\delta f_\varphi)^2}{C_0} & \dfrac{\delta f_\lambda+\delta^2f_r f_\varphi}{C_0}\\[3mm] \dfrac{\delta f_\varphi+\delta^2f_\lambda f_r}{C_0} & -\dfrac{\delta f_\lambda-\delta^2f_\varphi f_r}{C_0} & \dfrac{1+(\delta f_r)^2}{C_0} \end{array} \!\right), $$
and
\begin{eqnarray*} \alpha_{w1}&=&\dfrac{z_{\rm T}-z}{z_{\rm T}-z_{\rm s}}\phi_{\rm sx} ,\\ \alpha_{w2}&=&\dfrac{z_{\rm T}-z}{z_{\rm T}-z_{\rm s}}\phi_{\rm sy} ,\\ \alpha_{w3}&=&\dfrac{z_{\rm T}-z_{\rm s}}{z_{\rm T}} ,\\ \eta&=&1-\delta^2c_pZ_{\rm st}\dfrac{\partial\tilde{\Pi}}{\partial\hat{z}}\dfrac{\partial\tilde{\theta}}{\partial z} ,\\ C_0&=&(\delta f_\lambda)^2+(\delta f_\varphi)^2+\eta[1+(\delta f_r)^2] . \end{eqnarray*}
Notice that \(\hatw\) is the vertical velocity in the height-based terrain-following coordinate \(\hatz\). The relationship between w and \(\hatw\) is given by
\begin{equation} \label{eq28} w=\dfrac{z_{\rm T}-z}{z_{\rm T}-z_{\rm s}}\phi_{\rm sx}u+\dfrac{z_{\rm T}-z}{z_{\rm T}-z_{\rm s}}\phi_{\rm sy}v+\dfrac{z_{\rm T}-z}{z_{\rm T}}\hat{w} , \end{equation}
where z is the height level. Substitute w into the thermodynamic equation [Eq. (13)] to obtain θ' at the next time level:
\begin{equation} \label{eq29} \theta'_{n+1}=\xi_{\theta 1}\dfrac{\partial\Pi'}{\partial x}+\xi_{\theta 2}\dfrac{\partial\Pi'}{\partial y}+\xi_{\theta 3} \dfrac{\partial\Pi'}{\partial\hat{z}}+\xi_{\theta 0} , \end{equation}
where
\begin{eqnarray} \label{eq30} \xi_{\theta 1}&=&\dfrac{C_{31}\delta^2c_p\tilde{\theta}}{C_0}\dfrac{\partial\tilde{\theta}}{\partial z} ,\\ \label{eq31} \xi_{\theta 2}&=&\dfrac{C_{32}\delta^2c_p\tilde{\theta}}{C_0}\dfrac{\partial\tilde{\theta}}{\partial z} ,\\ \label{eq32} \xi_{\theta 3}&=&\dfrac{C_{31}Z_{\rm sx}+C_{32}Z_{\rm sy}+C_{33}Z_{\rm st}}{C_0}\delta^2c_p\tilde{\theta}\dfrac{\partial\tilde{\theta}}{\partial z} ,\\ \label{eq33} \xi_{\theta 0}&=&-\Bigg[\dfrac{C_{31}A_u+C_{32}A_v+C_{33}A_w}{C_0}-\nonumber\\[1mm] &&\dfrac{C_{33}\delta\left(s+c_pZ_{\rm st}\frac{\partial\tilde\Pi}{\partial\hat{z}}A_\theta\right)}{C_0}\Bigg] \delta\dfrac{\partial\tilde{\theta}}{\partial z}+A_\theta . \end{eqnarray}
When u,v,w and \(L_\Pi\) in Eq. (7) are substituted with Eqs. (14)-(16) and (12), a Helmholtz equation is deduced as
\begin{eqnarray} \label{eq34} \Pi'_{n+1}&\!=\!&\xi_{\Pi 1}\left[\left(\xi_{u1}\dfrac{\partial}{\partial x}+\xi_{u2}\dfrac{\partial}{\partial y}+ \xi_{u3}\dfrac{\partial}{\partial\hat{z}}\right)\Pi'+\xi_{u0}\right]+\nonumber\\ &&\xi_{\Pi 2}\left[\left(\xi_{v1}\dfrac{\partial}{\partial x}+\xi_{v2}\dfrac{\partial}{\partial y}+ \xi_{v3}\dfrac{\partial}{\partial\hat{z}}\right)\Pi'+\xi_{v0}\right]+\nonumber\\ &&\xi_{\Pi 3}\left[\left(\xi_{\hat{w}1}\dfrac{\partial}{\partial x}+\xi_{\hat{w}2}\dfrac{\partial}{\partial y}+\xi_{\hat{w}3} \dfrac{\partial}{\partial\hat{z}}\right)\Pi'+\xi_{\hat{w}0}\right]+\nonumber\\ &&\xi_{\Pi 4}\left\{\dfrac{\partial\left[\left(\xi_{u1}\frac{\partial}{\partial x}+\xi_{u2}\frac{\partial}{\partial y}+ \xi_{u3}\frac{\partial}{\partial\hat{z}}\right)\Pi'+\xi_{u0}\right]}{\partial x}\right.+\nonumber\\ &&\dfrac{\partial\left[\left(\xi_{v1}\frac{\partial}{\partial x}+\xi_{v2}\frac{\partial}{\partial y}+ \xi_{v3}\frac{\partial}{\partial\hat{z}}\right)\Pi'+\xi_{v0}\right]}{\partial y}+\nonumber\\ &&\left.\dfrac{\partial\left[\left(\xi_{\hat{w}1}\frac{\partial}{\partial x}\!+\!\xi_{\hat{w}2}\frac{\partial}{\partial y}\!+\! \xi_{\hat{w}3}\frac{\partial}{\partial\hat{z}}\right)\Pi'\!+\!\xi_{\hat{w}0}\right]}{\partial\hat{z}}\right\}\!+\!A_\Pi,\qquad \end{eqnarray}
where
\begin{eqnarray} \label{eq35} \xi_{\Pi 1}&=&-\delta\phi_{\rm sx}\left(\dfrac{z_{\rm T}-\hat{z}}{z_{\rm T}-z_{\rm s}}\dfrac{\partial\tilde{\Pi}}{\partial\hat{z}}- \dfrac{\tilde{\Pi}}{(\gamma-1)(z_{\rm T}-z_{\rm s})}\right),\\ \label{eq36} \xi_{\Pi 2}&=&-\delta\phi_{\rm sy}\left(\dfrac{z_{\rm T}-\hat{z}}{z_{\rm T}-z_{\rm s}}\dfrac{\partial\tilde{\Pi}}{\partial\hat{z}}- \dfrac{\tilde{\Pi}}{(\gamma-1)(z_{\rm T}-z_{\rm s})}\right),\\ \label{eq37} \xi_{\Pi 3}&=&-\delta\dfrac{\partial\tilde{\Pi}}{\partial\hat{z}} ,\\ \label{eq38} \xi_{\Pi 4}&=&-\dfrac{\delta\tilde{\Pi}}{\gamma-1} . \end{eqnarray}
The related 19 \(\Pi\) grid points, which serves the numerical solution of Eq. (35), are displayed in Fig. 1. It is worth noting that Eqs. (17)-(39) are all different from those in the original GRAPES model due to the full consideration of the 3D Coriolis term. The Helmholtz equation on the Yin-Yang grid is solved by using the GCR method with an incomplete LU factorization (ILU) (Liu and Cao, 2008) to speed up the convergence of the iterative algorithm. The classical Schwarz method is adopted to update the interface conditions of subdomains (Qaddouri et al., 2008). Bi-cubic Lagrangian interpolation is used for boundary updating. The threshold of absolute error, which is a sum of the Yin and Yang grids, is set to be 10-15 for the numerical convergence measurement.
3.2.1 Semi-Lagrangian transport on the Yin-Yang grid
Both the nonlinear terms and the departure-point-related terms are included in Ax in Eqs. (4)-(8). The departure point is calculated according to (Ritchie and Beaudoin, 1994). Halo regions are defined for each component zone to avoid multi-time data exchange during the parallel computation. Necessary data exchange is performed once per time step. When a departure point is located out of the computational region, it should be interpolated according to quantities in the halo zone (Fig. 2). The details are summarized as follows:
(1) Compute the position of the midpoint (Λ m,φ m,r m) at the half-time level on the sphere considering
\begin{eqnarray} \label{eq39} \lambda_{\rm m}&=&\lambda_{\rm a}-\dfrac{u_{\rm m}\Delta t}{2r_{\rm a}\cos(\varphi_{\rm a})} \left[1+\dfrac{\Delta t^2}{24r_{\rm a}^2}(u_{\rm m}^2\tan^2\varphi_{\rm a}-v_{\rm m}^2)\right],\quad\\ \label{eq40} \varphi_{\rm m}&=&\varphi_{\rm a}-\dfrac{v_{\rm m}\Delta t}{2r_{\rm a}}+\left(\dfrac{u_{\rm m}\Delta t}{2r_{\rm a}}\right)^2 \dfrac{\tan \varphi_{\rm a}}{2} ,\\ \label{eq41} r_{\rm m}&=&r_{\rm a}-\dfrac{w_{\rm m}\Delta t}{2} . \end{eqnarray}
where (Λ a,φ a,r a) represents the arrival point.
(2) Determine the velocity components at the midpoint (u m,v m,w m) with linear interpolation. If the departure point is located outside the computational domain, grid points in the halo region help accomplish the interpolation.
(3) Iterate (twice in this paper) the above two steps to modify the midpoint determination. The departure point (Λ d, φ d,r d) is then defined as
\begin{eqnarray} \label{eq42} \lambda_{\rm d}&=&\lambda_{\rm a}-\dfrac{u_{\rm m}\Delta t}{r_a\cos(\varphi_{\rm a})}\left(1-\dfrac{v_{\rm m}\Delta t}{2r_{\rm a}}\tan\varphi_{\rm a}\right),\\ \label{eq43} \varphi_{\rm d}&=&\varphi_{\rm a}-\dfrac{v_{\rm m}\Delta t}{r_{\rm a}}+\left(\sec^2\varphi_{\rm a}-\dfrac{2}{3}\right) \left(\dfrac{u_{\rm m}\Delta t}{2r_{\rm a}}\right)^2\dfrac{v_{\rm m}\Delta t}{2r_{\rm a}} ,\quad\\ \label{eq44} r_{\rm d}&=&r_{\rm a}-w_m\Delta t . \end{eqnarray}
Note that steps (2) and (3) are the same as in the original model, while step (3) is modified because of the existence of inner boundaries. Although the departure points can be out of the computational domain, they must be limited within the outer halo region boundaries. In this dynamic frame, four groups of departure points are calculated at scalar and vector grid points, separately, to define the coefficients in Eqs. (24), (28), (32), and (38).
3.2.2 Three-dimensional SISL integration for vectors on the Yin-Yang grid
The 3D SISL integration scheme (Qian et al., 1998) for vectors is used to calculate Au,v,w on the Yin-Yang grid
For scalars, the term \(A_\theta,\Pi\) can be interpolated at the departure point (denoted with superscript "*") directly:
\begin{equation} \label{eq45} A_x=(x')_\ast+[\alpha_\varepsilon\tilde{N}_x+\beta_\varepsilon(L_x+N_x)_\ast]\Delta t ,\quad x=\theta,\Pi , \end{equation}
where Lx(\(x=\theta,\Pi\)) is the linear term, Nx(\(x=\theta,\Pi\)) is the nonlinear variation and
\begin{eqnarray*} \beta_\varepsilon&=&1-\alpha_\varepsilon ,\\ \tilde{N}_{\rm x}&=&2(N_{\rm x})_n-(N_{\rm x})_{n-1} . \end{eqnarray*}
We note that the terms \(\xi_u0,v0,w0\) all contain Au,v,w,θ in Eqs. (20), (24) and (28). The prediction of 3D momentum calls for a good description of Au,v,w,θ, which depends on the computation of
\begin{equation} \label{eq46} Y_x=[x+\alpha_\varepsilon(L_x+N_x)\Delta t] ,\quad x=u,v,w,\theta , \end{equation}
at their departure points, denote by Yx*.
In the non-uniform vertical direction, a cubic Lagrangian interpolation,
\begin{equation} \label{eq47} Y_x^\ast(\hat{z})=\sum_{i=1}^4\dfrac{\Pi_{j=1,j\ne i}^4(\hat{z}_j-\hat{z})}{\Pi_{j=1,j\ne i}^4(\hat{z}_j-\hat{z}_i)}Y_x(\hat{z}_i) ,\quad x=u,v,w,\theta , \end{equation}
is used to ensure high-order accuracy. A second-order formulation for the vertical gradient of \(\Pi\) is given as
\begin{eqnarray} \label{eq48} \left.\dfrac{\partial\Pi}{\partial\hat{z}}\right|_{\hat{z}_{\Pi k}}&=&\Pi_{k-1}\dfrac{\hat{z}_k-\hat{z}_{k+1}}{(\hat{z}_{k-1}-\hat{z}_k)(\hat{z}_{k-1}-\hat{z}_{k+1})}+\qquad\nonumber\\ &&\Pi_k\dfrac{2\hat{z}_k-(\hat{z}_{k-1}+\hat{z}_{k+1})}{(\hat{z}_k-\hat{z}_{k-1})(\hat{z}_k-\hat{z}_{k+1})}+\nonumber\\ &&\Pi_{k+1}\dfrac{\hat{z}_k-\hat{z}_{k-1}}{(\hat{z}_{k+1}-\hat{z}_{k-1})(\hat{z}_{k+1}-\hat{z}_k)} . \end{eqnarray}
This ensures second-order accuracy without guaranteeing quantity conservation in the vertical transport. An alternative choice that achieves exact conservation is
\begin{equation} \label{eq49} \left.\dfrac{\partial\Pi}{\partial\hat{z}}\right|_{\hat{z}_{\Pi k}}=\dfrac{\Pi_{k+1}-\Pi_k}{(\hat{z}_{k+2}-\hat{z}_k)}+ \dfrac{\Pi_k-\Pi_{k-1}}{(\hat{z}_{k+1}-\hat{z}_{k-1})} , \end{equation}
but the accuracy decreases for non-uniform grids, where the \(\Pi_k\) is located at the middle level of \(\hatz_k\) and \(\hatz_k+1\). In the horizontal directions, four grids are needed in the halo region for a cubic Lagrangian interpolation at the departure point.