Advanced Search
Article Contents

An Improved Dynamic Core for a Non-hydrostatic Model System on the Yin-Yang Grid


doi: 10.1007/s00376-014-4120-5

  • A 3D dynamic core of the non-hydrostatic model GRAPES (Global/Regional Assimilation and Prediction System) is developed on the Yin-Yang grid to address the polar problem and to enhance the computational efficiency. Three-dimensional Coriolis forcing is introduced to the new core, and full representation of the Coriolis forcing makes it straightforward to share code between the Yin and Yang subdomains. Similar to that in the original GRAPES model, a semi-implicit semi-Lagrangian scheme is adopted for temporal integration and advection with additional arrangement for cross-boundary transport. Under a non-centered second-order temporal and spatial discretization, the dry nonhydrostatic frame is summarized as the solution of an elliptical problem. The resulting Helmholtz equation is solved with the Generalized Conjugate Residual solver in cooperation with the classic Schwarz method. Even though the coefficients of the equation are quite different from those in the original model, the computational procedure of the new core is just the same. The bi-cubic Lagrangian interpolation serves to provide Dirichlet-type boundary conditions with data transfer between the subdomains. The dry core is evaluated with several benchmark test cases, and all the tests display reasonable numerical stability and computing performance. Persistency of the balanced flow and development of both the mountain-induced Rossby wave and Rossby-Haurwitz wave confirms the appropriate installation of the 3D Coriolis terms in the semi-implicit semi-Lagrangian dynamic core on the Yin-Yang grid.
  • 加载中
  • Arakawa A., V. R. Lamb, 1997: Computational design of the basic dynamical processes of the UCLA General Circulation Model. Methods in Computational Physics, Vol. 17, J. Chang, Ed., Academic Press, San Francisco, USA, 173- 265.
    Baba Y., K. Takahashi, and T. Sugimura, 2010: Dynamical core of an atmospheric general circulation model on a Yin-Yang grid. Mon. Wea. Rev., 138, 3988- 4005.
    Charney J. G., N. A. Philips, 1953: Numerical integration of the quasi-geostrophic equations for barotropic and simple baroclinic Flow. J. Meteor., 10, 71- 99.
    Chen D., Coauthors, 2008: New generation of multi-scale NWP system (GRAPES): General scientific design. Chinese Sci. Bull., 53( 18), 3433- 3445.
    Jablonowski C., P. Lauritzen, R. Nair, and M. Taylor, 2008: Idealized test cases for the dynamical cores of atmospheric general circulation models: A proposal for the NCAR ASP 2008 summer colloquium. [Available online at http://esse.engin.umich. edu/admg/publications.php.]
    Kageyama A., T. Sato, 2004: "Yin-Yang grid": An overset grid in spherical geometry. Geochem. Geophys. Geosyst., 5, Q09005, doi: 10.1029/2004GC000734.
    Li X. L., D. H. Chen, X. D. Peng, F. Xiao, and X. S. Chen, 2006: Implementation of the semi-Lagrangian advection scheme on a quasi-uniform overset grid on a sphere. Adv. Atmos. Sci., 23(5), 792-801, doi: 10.1007/s00376-006-0792-9.
    Li X. L., D. H. Chen, X. D. Peng, K. Takahashi, and F. Xiao, 2008: A multimoment finite-volume shallow-water model on the Yin-Yang overset spherical grid. Mon. Wea. Rev., 136, 3066- 3086.
    Liu Y., J. W. Cao, 2008: ILU preconditioner for NWP system: GRAPES. Computer Engineering and Design, 29( 3), 731- 734. (in Chinese)
    McDonald A., J. R. Bates, 1989: Semi-Lagrangian integration of a grid point shallow water model on the sphere. Mon. Wea. Rev., 117, 130- 137.
    Peng X., F. Xiao, and K. Takahashi, 2006: Conservative constraint for a quasi-uniform overset grid on the sphere. Quart. J. Roy. Meteor. Soc., 132, 979- 996.
    Qian J., F. Semazzi, and J. S. Scroggs, 1998: A global nonhydrastatic semi-Lagrangian atmospheric model with orography. Mon. Wea. Rev., 126, 747- 771.
    Qaddouri A., L. Laayouni, J. Côté, and M. Gander, 2008: Optimized Schwarz methods with an overset grid for shallow-water equations: Preliminary results. Appl. Numer. Math., 58( 4), 459- 471.
    Ritchie H., C. Beaudoin, 1994: Approximations and sensitivity experiments with a baroclinic semi-Lagrangian spectral model. Mon. Wea. Rev., 122, 2391- 2399.
    Sadourny R. A., 1972: Conservative finite-difference approximations of the primitive equations on quasi-uniform spherical grids. Mon. Wea. Rev., 100, 136- 144.
    Sadourny R. A., A. Arakawa, and Y. Mintz, 1968: Integration of the nondivergent barotropic vorticity equation with an icosahedral-hexagonal grid for the sphere. Mon. Wea. Rev., 96, 351- 356.
    Satoh M., T. Inoue, and H. Miura, 2010: Evaluations of cloud properties of global and local cloud system resolving models using CALIPSO and CloudSat simulators. J. Geophys. Res., 115, D00H14, doi: 10.1029/2009JD012247.
    Skamarock W., J. Klemp, M. Duda, L. Fowler, and S. H. Park, 2012: A multi-scale nonhydrostatic atmospheric model using centroidal Voronoi tesselations and C-grid staggering. Mon. Wea. Rev., 140, 3090- 3105.
    Semazzi F., J. H. Qian, and J. S. Scroggs, 1995: A global nonhydrostatic, semi-Lagrangian, atmospheric model without orography. Mon. Wea. Rev., 123, 2534- 2550.
    Tomita H., M. Satoh, 2004: A new dynamical framework of nonhydrostatic global model using the icosahedral grid. Fluid Dyn. Res., 34, 357- 400.
    Williamson D. L., 2007: The evolution of dynamical cores for global atmospheric models. J. Meteor. Soc. Japan, 85B, 241- 269.
    Williamson D. L., J. B. Drake, J. J. Hack, R. Jackob, and P. N. Swarztrauber, 1992: A standard test for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102, 211- 224.
    Xue J. S., D. H. Chen, 2008: Design of GRAPES dynamical frame and the experiments. Scientific Design and Application of Numerical Prediction System GRAPES, J. H. Wang, Ed., Science Press, Beijing, 65- 136. (in Chinese)
  • [1] Qingchang QIN, Xueshun SHEN, Chungang CHEN, Feng XIAO, Yongjiu DAI, Xingliang LI, 2019: A 3D Nonhydrostatic Compressible Atmospheric Dynamic Core by Multi-moment Constrained Finite Volume Method, ADVANCES IN ATMOSPHERIC SCIENCES, 36, 1129-1142.  doi: 10.1007/s00376-019-9002-4
    [2] LI Xingliang, SHEN Xueshun, PENG Xindong, XIAO Feng, ZHUANG Zhaorong, CHEN Chungang, 2013: An Accurate Multimoment Constrained Finite Volume Transport Model on Yin-Yang Grids, ADVANCES IN ATMOSPHERIC SCIENCES, 30, 1320-1330.  doi: 10.1007/s00376-013-2217-x
    [3] LI Xingliang, CHEN Dehui, PENG Xindong, XIAO Feng, CHEN Xiongshan, 2006: Implementation of the Semi-Lagrangian Advection Scheme on a Quasi-Uniform Overset Grid on a Sphere, ADVANCES IN ATMOSPHERIC SCIENCES, 23, 792-801.  doi: 10.1007/s00376-006-0792-9
    [4] ZHAO Bin, ZHONG Qing, 2010: The Development of a Nonhydrostatic Global Spectral Model, ADVANCES IN ATMOSPHERIC SCIENCES, 27, 676-684.  doi: 10.1007/s00376-009-9080-9
    [5] PENG Xindong, CHANG Yan, LI Xingliang, XIAO Feng, 2010: Application of the Characteristic CIP Method to a Shallow Water Model on the Sphere, ADVANCES IN ATMOSPHERIC SCIENCES, 27, 728-740.  doi: 10.1007/s00376-009-9148-6
    [6] Shi Yong, Jiang Weimei, 1998: The Numerical Simulation on the PBL Structure and Its Evolution over Small-Scale Concave Terrain, ADVANCES IN ATMOSPHERIC SCIENCES, 15, 99-106.  doi: 10.1007/s00376-998-0021-9
    [7] Wang Yunfeng, Wu Rongsheng, Wang Yuan, Pan Yinong, 1999: Application of Variational Algorithms in Semi-Lagrangian Framework, ADVANCES IN ATMOSPHERIC SCIENCES, 16, 419-430.  doi: 10.1007/s00376-999-0020-5
    [8] Pei HUANG, Chungang CHEN, Xingliang LI, Xueshun SHEN, Feng XIAO, 2022: An Adaptive Nonhydrostatic Atmospheric Dynamical Core Using a Multi-Moment Constrained Finite Volume Method, ADVANCES IN ATMOSPHERIC SCIENCES, 39, 487-501.  doi: 10.1007/s00376-021-1185-9
    [9] Yang XIA, Bin WANG, Lijuan LI, Li LIU, Jianghao LI, Li DONG, Shiming XU, Yiyuan LI, Wenwen XIA, Wenyu HUANG, Juanjuan LIU, Yong WANG, Hongbo LIU, Ye PU, Yujun HE, Kun XIA, 2024: A Neural-Network-Based Alternative Scheme to Include Nonhydrostatic Processes in an Atmospheric Dynamical Core, ADVANCES IN ATMOSPHERIC SCIENCES.  doi: 10.1007/s00376-023-3119-1
    [10] WANG Xiaocong, LIU Yimin, WU Guoxiong, Shian-Jiann LIN, BAO Qing, 2013: The Application of Flux-Form Semi-Lagrangian Transport Scheme in a Spectral Atmosphere Model, ADVANCES IN ATMOSPHERIC SCIENCES, 30, 89-100.  doi: 10.1007/s00376-012-2039-2
    [11] Jie TANG, Chungang CHEN, Xueshun SHEN, Feng XIAO, Xingliang LI, 2021: A Positivity-preserving Conservative Semi-Lagrangian Multi-moment Global Transport Model on the Cubed Sphere, ADVANCES IN ATMOSPHERIC SCIENCES, 38, 1460-1473.  doi: 10.1007/s00376-021-0393-7
    [12] HUANG Bo, CHEN Dehui, LI Xingliang, LI Chao, , 2014: Improvement of the Semi-Lagrangian Advection Scheme in the GRAPES Model: Theoretical Analysis and Idealized Tests, ADVANCES IN ATMOSPHERIC SCIENCES, 31, 693-704.  doi: 10.1007/s00376-013-3086-z
    [13] Lucas HARRIS, 2021: A New Semi-Lagrangian Finite Volume Advection Scheme Combines the Best of Both Worlds, ADVANCES IN ATMOSPHERIC SCIENCES, 38, 1608-1609.  doi: 10.1007/s00376-021-1181-0
    [14] Chen Jiabin, Wang Jun, 1996: Studies on Non-interpolating Semi-Lagrangian Scheme and Numerical Solution to KdV Equation, ADVANCES IN ATMOSPHERIC SCIENCES, 13, 265-268.  doi: 10.1007/BF02656869
    [15] YUAN Zhuojian, JIAN Maoqiu, 2003: A Linear Diagnostic Equation for the Nonhydrostatic Vertical Motion W in Severe Storms, ADVANCES IN ATMOSPHERIC SCIENCES, 20, 875-881.  doi: 10.1007/BF02915511
    [16] SHEN Xueshun, Akimasa SUMI, 2005: A High Resolution Nonhydrostatic Tropical Atmospheric Model and Its Performance, ADVANCES IN ATMOSPHERIC SCIENCES, 22, 30-38.  doi: 10.1007/BF02930867
    [17] Yufan DAI, Qingqing LI, Xinhang LIU, Lijuan WANG, 2023: A Lagrangian Trajectory Analysis of Azimuthally Asymmetric Equivalent Potential Temperature in the Outer Core of Sheared Tropical Cyclones, ADVANCES IN ATMOSPHERIC SCIENCES, 40, 1689-1706.  doi: 10.1007/s00376-023-2245-0
    [18] Liu Yudi, Ji Zhongzhen, Wang Bin, 2002: Study on Computational Properties of Several Vertical Grids with a Nonhydrostatic Model in Comparison to Analytical Solutions, ADVANCES IN ATMOSPHERIC SCIENCES, 19, 528-543.  doi: 10.1007/s00376-002-0084-y
    [19] Sun Litan, Huang Meiyuan, 1994: Improving the Vorticity-Streamfunction Method to Solve Two-Dimensional Anelastic and Nonhydrostatic Model, ADVANCES IN ATMOSPHERIC SCIENCES, 11, 247-249.  doi: 10.1007/BF02666551
    [20] Syed Faizan Haider, K. H. Asif, Amjad Hussain Gilani, 1992: Weather Yield Model for the Semi Tropical Region (Pakistan), ADVANCES IN ATMOSPHERIC SCIENCES, 9, 367-372.  doi: 10.1007/BF02656947

Get Citation+

Export:  

Share Article

Manuscript History

Manuscript received: 11 June 2014
Manuscript revised: 22 September 2014
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

An Improved Dynamic Core for a Non-hydrostatic Model System on the Yin-Yang Grid

  • 1. State Key Laboratory of Severe Weather, Chinese Academy of Meteorological Sciences, Beijing 100081
  • 2. Center of Numerical Weather Prediction, China Meteorological Administration, Beijing 100081

Abstract: A 3D dynamic core of the non-hydrostatic model GRAPES (Global/Regional Assimilation and Prediction System) is developed on the Yin-Yang grid to address the polar problem and to enhance the computational efficiency. Three-dimensional Coriolis forcing is introduced to the new core, and full representation of the Coriolis forcing makes it straightforward to share code between the Yin and Yang subdomains. Similar to that in the original GRAPES model, a semi-implicit semi-Lagrangian scheme is adopted for temporal integration and advection with additional arrangement for cross-boundary transport. Under a non-centered second-order temporal and spatial discretization, the dry nonhydrostatic frame is summarized as the solution of an elliptical problem. The resulting Helmholtz equation is solved with the Generalized Conjugate Residual solver in cooperation with the classic Schwarz method. Even though the coefficients of the equation are quite different from those in the original model, the computational procedure of the new core is just the same. The bi-cubic Lagrangian interpolation serves to provide Dirichlet-type boundary conditions with data transfer between the subdomains. The dry core is evaluated with several benchmark test cases, and all the tests display reasonable numerical stability and computing performance. Persistency of the balanced flow and development of both the mountain-induced Rossby wave and Rossby-Haurwitz wave confirms the appropriate installation of the 3D Coriolis terms in the semi-implicit semi-Lagrangian dynamic core on the Yin-Yang grid.

1. Introduction
  • The dynamic core is the heart of the atmospheric model. It determines the computing characteristics such as the numerical accuracy and computational efficiency. Use of a high-order spatial finite differencing scheme, semi-Lagrangian transport, advanced temporal integration, and other state-of-the-art techniques has improved model representation significantly. In addition, with the rapid development of computer systems, there is increasing demand for global high-resolution numerical weather predictions. With a conventional latitude-longitude grid, the difference between the mesh size at the equator and in the polar regions becomes larger as the model resolution increases. The large and complex computations involved in running global high-resolution models calls for a quasi-uniform grid on a sphere. For example, the Nonhydrostatic Icosahedral Atmospheric Model (NICAM, Tomita and Satoh, 2004) and Model for Prediction Across Scales-Atmosphere (MPAS-A, Skamarock et al., 2012) have been developed recently for multi-scale simula- tions and cloud-resolving forecasts. Mesoscale phenomena and even cloud development can be resolved with a global atmospheric model (Satoh et al., 2010). A global non-hydrostatic model system, the Global/Regional Assimilation and Prediction System (GRAPES, Xue and Chen, 2008), was developed on the latitude-longitude grid system at the Chinese Meteorological Administration. It is a fully compressible model designed for numerical weather prediction across scales. The height-based terrain-following coordinate is used to incorporate orography for the convenience of bottom boundary arrangement. The GRAPES dynamic core is solved with the semi-implicit semi-Lagrangian (SISL) method containing the spherical departure-point determination (Ritchie and Beaudoin, 1994) and the non-centered finite differencing scheme (Semazzi et al., 1995). The 3D vector SISL scheme (Qian et al., 1998) is adopted to discretize the momentum equations. The Charney-Phillips grid (Charney and Philips, 1953) and Arakawa-C grid (Arakawa and Lamb, 1997) were chosen for the staggering of physical variables in the vertical and horizontal directions, respectively.

    In high-resolution cases, however, numerical approaches developed for the longitude-latitude coordinate face additional difficulties, such as pole singularity and the convergence of meridians in the polar regions. The meridians converge closer to the poles, which makes the latitudinal grid interval null at the pole points. The singularity displays a problem in vector representation. On the other hand, the singularity is not a physical property, but a problem of coordinate system selection. The problem can be overcome in a variety of ways, such as diagnosing the polar wind with the minimization principle (McDonald and Bates, 1989). In a numerical model, the time step is rigidly limited by the smallest grid spacing at the poles. This shortcoming turns into a serious issue for effective and economic model integrations in high-resolution models. In fact, the problem becomes worse and more complicated in the case of the SISL model. The existence of tan φ and sec φ in the horizontal advection terms reduces the accuracy of computation in high latitude areas when dealing with the spherical departure-point in the Ritchie scheme. In a latitude-longitude mesh model, subgrid-scale processes at midlatitudes may be resolved as grid-scale processes in the polar regions due to the relatively high resolution there, which confuses the physical interaction among scales. Different physics schemes are then required in the corresponding model for proper simulations.

    Owing to the aforementioned problems, model development using quasi-uniform grids on a sphere is now an important topic for multiscale simulation and numerical weather forecasting. The most familiar designs of quasi-uniform grids are the icosahedral grid (Sadourny et al., 1968), the cubed grid (Sadourny, 1972), and the Yin-Yang grid (Kageyama and Sato, 2004). The advantages and disadvantages of each of these grids have been discussed in (Williamson, 2007), in which the Yin-Yang grid was selected as an appropriate successor to the longitude-latitude coordinate. Most numerical algorithms based on the longitude-latitude coordinate can be straightforwardly used on the Yin-Yang grid without any change. Mesh refinement and the practice of grid nesting are convenient because of the structured and regular grid distribution. Global and regional models can share the same dynamic frame in the Yin-Yang grid system. The disadvantage of the Yin-Yang grid, however, is the existence of inner boundaries. In this paper, we show the details of the computing procedure of the SISL method when a 3D Coriolis forcing is added. Redefinition of the coefficients of the Helmholtz equation and arrangement of the boundary consideration in the Generalized Conjugate Residual (GCR) iteration are displayed with respect to the dynamic core on the Yin-Yang grid.

    The paper is organized as follows: A brief description of the Yin-Yang grid with the interpolation scheme for inner boundaries is presented in next section. The detailed computing procedure of the nonhydrostatic dynamic core of GRAPES on the Yin-Yang grid (GRAPES_YY) follows in section 3, including the SISL method, computing of the corresponding coefficients, and arrangement of cross-boundary transport. Numerical results of some standard tests are presented in section 4, and a conclusion to the study is provided in section 5.

2. The Yin-Yang grid
  • The Yin-Yang grid, first proposed by (Kageyama and Sato, 2004), is composed of two identical component zones. The pair of zones are combined in a complementary way to cover the sphere with overlaps at their boundaries. Each component zone is a low latitude part of the longitude-latitude grid, and one is rotated by 90° to fit with the other. The grid spacing is quasi-uniform——the minimum/maximum ratio of the grid spacing is about 0.707——and the coordinate is orthogonal without any singularity. Therefore, high-order interpolation schemes and finite difference methods that have been developed on the longitude-latitude coordinate can be applied on the Yin-Yang grid. In addition, a set of parallel computation methods can be easily introduced into the system thanks to the identical structure of the two zones. Like any other overset grid, however, the Yin-Yang grid requires interpolation at the boundaries, which might reduce the accuracy of numerical integrations and frequent communications between processors will definitely decrease the computing speed of the parallel program. Global conservation is another issue for the overlapped Yin-Yang grid. The existence of inner boundaries becomes a problem for the mass-flux balance between the component zones. (Peng et al., 2006) developed a numerical constraint to guarantee that the fluxes at the boundaries of the two components are identical, which achieves local and global conservation.

    To solve the boundary problem of the Yin-Yang grid, a 2D Lagrange interpolation was introduced for boundary data exchange in Li et al. (2006, 2008) and (Baba et al., 2010). In their work, the results of benchmark tests showed that the presence of the overset region does not significantly affect the dynamics on both long and short time scales when the high-order interpolation methods are applied. In a Yin-Yang grid system, boundary data exchange occurs in the halo region that is aligned to the bounds of each domain. Quantities in the halo region are not updated with temporal integration of prognostic equations, but with interpolation from another zone. In this paper, two overset grids are defined for the cubic Lagrange interpolation. The coordinate conversion between Yin and Yang grids can be expressed as

    \begin{equation} \label{eq1} \left\{ \begin{array}{l} \cos\varphi_o\cos\lambda_o=-\cos\varphi_e\cos\lambda_e\\[1.5mm] \cos\varphi_o\sin\lambda_o=\sin\varphi_e\\[1.5mm] \sin\varphi_o=\cos\varphi_e\sin\lambda_e \end{array} \right., \end{equation}

    where Λ and φ represent longitude and latitude, respectively, and the subscripts e and o denote the Yin and Yang grids. The scalar quantity can be interpolated directly with bi-cubic Lagrangian interpolation (Li et al., 2006), and vector transformation from the Yang to Yin zone gives

    \begin{equation} \label{eq2} \left( \begin{array}{c} v_e \\[1.5mm] u_e \end{array} \right)=\left( \begin{array}{c@{\quad}c} -\sin\lambda_o\sin\lambda_e & -\cos\lambda_o/\cos\varphi_e\\[1.5mm] \cos\lambda_o/\cos\varphi_e & -\sin\lambda_o\sin\lambda_e \end{array} \right)\left( \begin{array}{c} v_o\\[1.5mm] u_o \end{array} \right), \end{equation}

    and vice versa.

3. Dynamical frame of the GRAPES_YY
  • 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.

    Figure 1.  Distribution of the 19 \(\Pi\) grid points involved in the discretization of the Helmholtz equation.

  • 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.

    Figure 2.  Determination of the departure point and calculation of Y* on the Arakawa-C staggered grid.

    (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.

4. Numerical results of benchmark tests
  • For the validation of the dynamical modification concerning 3D Coriolis and trajectory computation across boundaries, several benchmark tests are carried out to check the computational accuracy and the performance of GRAPES_ YY. The model top is defined as 32.5 km, and the model atmosphere is divided into 36 non-uniform levels. The horizontal resolution is 2.5°. There is no viscosity added in the dry core.

  • This test is a 3D extension of Test 2 in (Williamson et al., 1992). The initial state is defined as

    \begin{eqnarray} \label{eq50} &&u=u_0(\cos\varphi\cos\alpha+\cos\lambda\sin\varphi\sin\alpha) ,\quad \\ \label{eq51} &&v=-u_0\sin\varphi\sin \alpha ,\\ \label{eq52} &&w=0 ,\\ \label{eq53} &&c_p\theta\dfrac{\partial\Pi}{\partial z}=-g ,\\ \label{eq54} &&\dfrac{c_p\theta}{r}\dfrac{\partial\Pi}{\partial\varphi}=-2\Omega u\sin\varphi-\dfrac{u^2}{r}\tan\varphi , \end{eqnarray}

    where u0 is set to be 20 m s-1 and α is flow orientation angle, which is 0 here. Thirty-day integration results of \(\Pi',u,v,w\) and the differences between the numerical solution and the exact one, with a time step of 1800 s, are shown in Figs. 3a-d. Corresponding results from GRAPES are presented in Figs. 3e-h for comparison. The \(\Pi'\) in Fig. 3a shows its zonal parallel contours after the 30-day integration, while the absolute error is about -1.0× 10-4 at the equator and 5.0×10-5 in midlatitude regions. Zonal wind (Fig. 3b) keeps its initial state with a perturbation of -0.12 m s-1 at the boundary of the Yin-Yang grid. Meridional and vertical winds, which display as absolute error in Figs. 3c and d, show their order of about 10-3 and 10-7 m s-1, respectively. Meanwhile, errors of meridional and vertical winds reach 0.1 and 2× 10-3 m s-1 in the polar regions of the GRAPES model (Figs. 3g and h). Relatively small error is found with the new dynamic core on the Yin-Yang grid. It is clear that obvious numerical error exists at the boundaries of the overset grid in GRAPES_YY, and the error appears in the pole regions in the original GRAPES core. The new core shows a much smaller error than the original one. Owing to fact the analytical solutions of the meridional and vertical velocities are null, numerical errors appear as highly significant. Consequently, a boundary trail is revealed in Figs. 3c and d. The time series of the corresponding error norms \(\ell_1,\ell_2\) and \(\ell_\infty\) of both the scalar and vectors are given for GRAPES_YY in Fig. 4. Error norms increase with time. The \(\ell_1\) and \(\ell_2\) norms of the \(\Pi'\) are 0.002 and 0.0025 at day 30, while those of the velocity are 0.009 and 0.0095, respectively. This numerical test confirms the stability of the SISL method and the proper installation of the 3D Coriolis force in the nonhydrostatic frame.

    Figure 3.  Numerical results of the steady-state zonal flow test at day 30: (a) perturbation of Exner function \(\Pi'\) (contours) and the error (shaded, 10-4), (b) u (contours, m s-1) and the error (shaded), (c) v (10-2 m s-1), and (d) vertical wind (10-7 m s-1) with the GRAPES_YY new core in comparison with (e) \(\Pi'\) (contours) and the error (shaded, 10-4), (f) u (contours, m s-1) and the error (shaded, m s-1), (g) v (10-1 m s-1), and (h) vertical wind (10-3 m s-1) with the original GRAPES model.

    In the zonal flow case, small interpolation error at the boundaries is clearly displayed with the meridional wind and vertical velocity. We find the error of v and w to be 10-3 and 10-7, respectively, which is negligible in comparison with the u component. But does the error destroy the model stability in a non-zero meridional wind case? We also show the numerical results of the balanced flow with an angle of 45° in Fig. 5. This configuration seems to be the harshest for the Yin-Yang grid because of the orthogonality of the two sub-zones. Owing to the zonal wind enhancement at high latitudes in this test case, the time step will be tightly limited in a latitude-longitude grid system. On the quasi-uniform Yin-Yang grid, however, the time step remains the same as in the former. A reasonable distribution is found with all the scalar and vector quantities. The scalar \(\Pi'\), zonal wind and vertical wind display equivalent computational error as in Fig. 3. Even though the meridional wind shows larger error because of its enhancement, the boundary trail is nearly invisible in the vector fields.

    Figure 4.  The same as Fig. 3 but for the flow with an orientation angle α of 45° with the GRAPES_YY core. The error in (c) is plotted in m s-1.

    Figure 5.  Error norms \(\ell_1\), \(\ell_2\) and \(\ell_\infty\) of (a) \(\Pi'\) and (b) velocity in the steady-state flow test during the 30-day integration.

  • In this test, the dynamic core with 3D Coriolis force and the SISL solver is tested with topography. The initial wind velocity is the same as the previous one. The mountain height is given by

    \begin{equation} \label{eq55} h=h_0\exp\left[-\left(\dfrac{D}{R}\right)^2\right] , \end{equation}

    where h0=2000 m determines the peak height of the mountain and R=1500 km is the mountain half width; D, the distance to the mountain center c c)=(π/2,π/6):

    \begin{equation} \label{eq56} D=a\arccos[\sin\varphi_c\sin\varphi+\cos\varphi_c\cos\varphi\cos(\lambda-\lambda_c)] , \end{equation}

    The pressure field is initially given a hydrostatic balance state,

    \begin{equation} \label{eq57} p_{\rm s}(\lambda,\varphi)\!=\!p_{\rm sp}\exp\left[-\dfrac{aN^2u_0}{2g^2\kappa}\left(\dfrac{u_0}{a}\!+\!2\Omega\right)(\sin^2\varphi-1)\!-\!\dfrac{N^2h}{g\kappa}\right], \end{equation}

    where N=0.0182 s-1 is the Brunt-Valsala frequency, \(\kappa=2/7,u_0=20\) m s-1, a=6371.229 km is mean radius of earth, Ω is earth's angular velocity and p sp=930 hPa denotes the surface pressure at the South Pole. No analytical solution is available in this test case, but the results with a spectral method (Jablonowski et al., 2008) can be referenced. Geopotential height, temperature, and the horizontal wind components u and v at the 700 hPa level of day 15 are given in Fig. 6 in comparison to the results of the original GRAPES. All the results of the new core integration are comparable with those in (Jablonowski et al., 2008), even though a low-resolution configuration is used here. The figures illustrate a proper evolution of the mountain-induced Rossby wave, and no boundary trail is found with the non-zero velocity in this case. The original model, however, displays the Rossby wave as not well developed, due to the low-resolution configuration. Therefore, we can again confirm good numerical performance via this test case.

    Figure 6.  Numerical result of the mountain-wave test at 700 hPa with the GRAPES_YY (contours) and original GRAPES (shaded) model: (a) geopotential height (gpm), (b) temperature (°C), (c) u (m s-1) and (d) v (m s-1) at day 15.

  • The initial velocity field is given by

    \begin{eqnarray} u&=&aM\cos\varphi\!+\!aK\cos^{n-1}\varphi\cos(c\lambda)(c\sin^2\varphi\!-\!\cos^2\varphi) ,\quad\ \ \\[0.5mm] \label{eq59} v&=&-aKc\cos^{n-1}\varphi\sin\varphi\sin(c\lambda) ,\\[0.5mm] \label{eq60} w&=&0, \end{eqnarray}

    where c=4 denotes the wave number, and M=K≈1.962× 10-6 s-1. Profiles of the temperature and pressure are given by

    \begin{eqnarray} \label{eq61} T&=&T_0-\Gamma\left(z-\dfrac{\widetilde{\Phi}}{g}\right),\\ \label{eq62} p&=&p_{\rm ref}\left(\dfrac{T}{T_0}\right)^{\frac{g}{r_{\rm d}\Gamma}} , \end{eqnarray}

    where p ref=955 hPa and T0=288 K with the moist-adiabatic lapse rate Γ=0.0065 km-1. The geopotential is given as

    \begin{equation} \label{eq63} \Phi(\lambda,\varphi)=a^2A(\varphi)+a^2B(\varphi)\cos(c\lambda)+a^2C(\varphi)\cos(2c\lambda), \end{equation}

    where

    \begin{eqnarray} \label{eq64} A(\varphi)&\!=\!&\dfrac{M(2\Omega+M)}{2}\cos^2\varphi+\dfrac{K^2}{4}\cos^{2n}\varphi[(c+1)\cos^2\varphi+\nonumber \end{eqnarray} \begin{eqnarray} &&(2c^2-c-2)]-\dfrac{c^2K^2}{2}\cos^{2(c-1)}\varphi ,\\[1.5mm] \label{eq65} B(\varphi)&\!=\!&\dfrac{2(\Omega+M)K}{(c\!+\!1)(c\!+\!2)}\cos^c\varphi[(c^2+2c+2)\!-\!(c+1)^2\cos^2\varphi] ,\nonumber\\\\[1mm] \label{eq66} C(\varphi)&\!=\!&\dfrac{K^2}{4}\cos^{2c}\varphi[(c+1)\cos^2\varphi-(c+2)] . \end{eqnarray}

    The time step is changed to 600 s in this test for the serious limitation of linear computational stability. The numerical results of the geopotential height, u, and v at 500 hPa and surface pressure at day 14 are plotted in Fig. 7. The four-wave structure propagates correctly in geopotential height, surface pressure and the horizontal wind field in this low-resolution model. No obvious numerical deformation of the wave is observed, and the wave displays smooth propagation at the Yin-Yang boundaries with the classic Schwarz scheme.

    Figure 7.  Numerical result of the Haurwitz-Rossby wave at 500 hPa with the new dynamic core on the Yin-Yang grid: (a) geopotential height (gpm), (b) surface pressure (hPa), (c) u (m s-1) and (d) v (m s-1) at day 14.

    The cost of the classic Schwarz solver is about 24.92% of the model total expense due to the frequent information exchange for the boundary constraint of the Helmholtz equation. Of course, the cost varies with the iteration in the GCR solver. With the help of the ILU preconditioner, the convergence of the GCR solver shows great efficiency. The iteration before its convergence is listed in Fig. 8 for the first 100-step integration. Rapid convergence of the solvers is achieved for the overset grid, and the iteration tends to decrease with time.

    Figure 8.  Iterations of the GCR and classic Schwarz solvers before convergence in the first 100 steps of the 45° steady-state flow, zonal flow over a mountain, and the Rossby-Haurwitz wave test cases.

    The test of the zonal flow over a mountain shows more iterations than the others for its strong time-dependent current. Conclusion An improved dynamic core of the GRAPES model is successfully reconstructed on a quasi-uniform Yin-Yang grid, which is free of pole singularity. Three-dimensional Coriolis force has been introduced to the new frame, which makes the code identical between the Yin and Yang components. The departure point across boundaries is fixed with the help of the halo region, and the SISL scheme for 3D vectors is implemented into the dynamical core on the Yin-Yang grid. Numerical results of 3D benchmark tests reveal strong computational stability and reasonable performance. The results also show the property of the 3D Coriolis installation and the reconstruction of the Helmholtz equation in the SISL integration. The new nonhydrostatic core displays reasonable numerical results in three idealized tests with or without topography. The classic Schwarz method, which updates the boundary with a bi-cubic Lagrangian interpolation, is generally efficient for the constraint of global convergence of the numerical solution. On the other hand, relatively expensive cost and numerical oscillations at the boundary are also observed in the SISL integration. Further investigation on the interpolation procedure of the determination of coefficients for the Helmholtz equation in the classical Schwarz method is demanded by the SISL integration of non-hydrostatic models on the Yin-Yang grid. Additional developments, such as the boundary constraint, parallelization for large computations, and consideration of vapor and the installation of physics, will be pursued in future work.

Reference

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return