Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1225-6692(Print)
ISSN : 2287-4518(Online)
Journal of the Korean earth science society Vol.39 No.4 pp.317-326

Effect of Nonuniform Vertical Grid on the Accuracy of Two-Dimensional Transport Model

Chung-Hui Lee1, Hyeong-Bin Cheong1*, Hyun-Ju Kim1, Hyun-Gyu Kang1,2
1Department of Environmental Atmospheric Sciences, Pukyong National University, Busan 48513, Korea
2BK21 plus Project of the Graduate School of Earth Environmental Hazard System, Busan 48513, Korea
Corresponding author: Tel: +82-51-629-6642 Fax: +82-51-629-6636
July 25, 2018 August 16, 2018 August 23, 2018


Effect of the nonuniform grid on the two-dimensional transport equation was investigated in terms of theoretical analysis and finite difference method (FDM). The nonuniform grid having a typical structure of the numerical weather forecast model was incorporated in the vertical direction, while the uniform grid was used in the zonal direction. The staggered and non-staggered grid were placed in the vertical and zonal direction, respectively. Time stepping was performed with the third-order Runge Kutta scheme. An error analysis of the spatial discretization on the nonuniform grid was carried out, which indicated that the combined effect of the nonuniform grid and advection velocity produced either numerical diffusion or numerical adverse-diffusion. An analytic function is used for the quantitative evaluation of the errors associated with the discretized transport equation. Numerical experiments with the non-uniformity of vertical grid were found to support the analysis.



    Numerical weather prediction (NWP) model uses various discretization method and grid structure (Haltiner and Williams, 1980; Gates, 1992; Klemp, 2011; Kang and Cheong, 2017). As the computing resources to run the model becomes powerful, the spatial resolution of the NWP model has been increased steadily (Tomita and Satoh, 2004). In general, the vertical resolution is much higher than the horizontal resolution because the vertical variation of atmospheric motion and state is much stronger than the horizontal due to a large horizontal-to-vertical ratio of atmospheric layer (Haltiner and Williams, 1980).

    In traditional NWP models the grid-size of the horizontal grid is not uniform, as is characterized by the latitude-longitude grid system (Smagorinsky, 1963; Arakawa, 1972; Hoskins and Simmons, 1975; Arakawa and Lamb, 1977). Although the lat-lon grid system is still used in some models and operational NWP centers (Saito et al., 2006; Skamarock and Klemp, 2008), the quasi-uniform grid system has gained increasing popularity (Rivier et al., 2002; Tomita and Satoh, 2004; Nair et al., 2005; Zhang and Rancic, 2007; Marras et al., 2015; Choi and Hong, 2016; Kang and Cheong, 2017, 2018). In the vertical direction, the grid-size varies substantially with height, typically decreasing from the middle level towards lower and upper boundaries (Gal-Chen and Somerville, 1975; Haltiner and Williams, 1980; Simmons and Burridge, 1981; Saito et al., 2006). This is because, unlike the horizontal domain, the atmosphere is bounded by the Earth’s surface and upper boundary, and it is needed to resolve the atmospheric motion and state with finer grid near the boundaries where the flow speed changes steeply (Haltiner and Williams, 1980; Leuenberger et al., 2010; Klemp, 2011). Since there is no clear physical boundary in the upper atmosphere, not as the surface boundary, the NWP model usually sets a constant pressure level or height level as the upper boundary.

    The nonuniform grid in the vertical direction is regarded as an essential component of the NWP model as is convinced in numerical weather predictions (Haltiner and Williams, 1980; Saito et al., 2006; Klemp 2011). However, the nonuniform grid, in comparison with the uniform grid, is known to provide less discretization accuracy for the governing equations of the atmospheric motions and state (Arakawa and Lamb, 1977; Durran, 1999). In spite of this inherent deficiency of the nonuniform grid, little attention has been paid to the systematic accuracy assessment so far.

    In this study, the accuracy of the nonuniform grid for the two-dimensional transport equation is examined. For simplicity and clarity, the grid structure is given to be symmetric with respect to the middle point of the vertical domain. The nonuniform grid is established in such a way to maintain the ratio of larger-to-smaller grid-size of neighboring vertical grid-points to be the same. The paper is organized as follows. Next section will present the two-dimensional transport equation and boundary conditions. Section 3 shows the discretization method and theoretical error analysis with a focus on the nonuniformity of the grid-size. Initial conditions and analytic solution are provided in section 4. Results of numerical experiment and discussion in comparison with the theoretical analysis will be given in the subsequent section. Concluding remark is presented in the final section.

    Transport Equation and Boundary Conditions

    Two dimensional transport (advection) equation for an arbitrary scalar variable h is written as:

    h t + u h x + w h z = 0 ,

    where u and w are the advection velocity of x and z direction, respectively. Velocity field is given to be nondivergent, that is, u x + w z = 0 . The nondivergence of the flow makes it feasible to introduce the streamfunction such that u = ψ z and w = ψ x .

    Model domain is set as 0 x L and 0 z H , and the boundary condition of Eq. (1) is set as follows:

    h ( 0 , z ) = h ( L , z ) u ( 0 , z ) = u ( L , z ) w ( 0 , z ) = w ( L , z ) w ( x , 0 ) = w ( x , H ) = 0 u z ( x , 0 ) = u z ( x , H ) = 0.

    In the case where the scalar variable at boundaries is given initially zero, it remains as zero-valued at any time during time evolution. This is expressed as the following relationship:

    h ( x , 0 ) = h ( x u t , 0 ) h ( x , H ) = h ( x u t , H ) .

    Discretization Method

    Equation (1) is solved numerically for a given initial field of h and velocity field (u, w) using the third-order Runge Kutta scheme (Nair et al., 2005) and the FDM for the time and the space, respectively. The grid system used in this study is uniform in xdirection but nonuniform in the vertical direction, as shown in Fig. 1. Such a grid system is a typical one adopted in the numerical weather prediction models (Haltiner and Williams, 1980; Tomita and Satoh, 2004; Skamarock and Klemp, 2008). The grid interval in z-direction decreases monotonically from the midlevel toward the boundaries. In the present study, the ratio of two adjacent grid-intervals is kept at a constant value γ in the lower half domain and γ −1 in the upper half domain, which is controlled as a model parameter. A typical value of γ in operational model is 1.0-1.2 (e.g., Saito et al., 2006; Haiden et al., 2017). Location of variables is depicted in Fig. 2 together with the vertical nonuniform grids. Non-staggered grid is used in the horizontal direction, while staggered grid is adopted in the vertical direction: All variables are located in the same vertical line, and the scalar variable and horizontal velocity are placed in the middle of vertical grid-points. By choosing this gridsystem, it is possible to focus on the effect of the vertical nonuniform grid on the transport equation.

    Third-order Runge Kutta scheme is written as follows (Nair et al., 2005):

    h 1 = h ( n ) + δ t T ( h ( n ) ) h 2 = 3 4 h ( n ) + 1 4 h 1 δ t 4 T ( h 1 ) h n + 1 = 1 3 h ( n ) + 3 4 h 2 δ t 4 T ( h 2 ) ,

    where the superscript n means time-level of timestepping procedure, that is, t n ( = n δ t ) with δt being the time-step size, and T implies the time tendency of the equation.

    Spatial discretization is performed with FDM: For the horizontal discretization the second-order centered FDM is used, while for the vertical discretization the staggered-grid 3-points stencil FDM is chosen. The discrete form of these methods can be derived based on the Taylor series. The procedure to derive the finite difference form related to these methods is explained briefly. Taylor series of an arbitrary function f ( z ) at z + α is given as follows

    f ( z + α ) = f ( z ) + f α + f α 2 2 ! + f ( 3 ) α 3 3 ! + ... + f ( n ) α n n ! + ... = n = 0 n = f ( n ) α n n ! ,

    where f ( n ) = n f z n and n ! = n ( n 1 ) 2 1 . Let the function f ( z ) be defined on discrete grid-points z k , which is denoted by f ( z k ) . The grid interval is defined as Δ k = z k z k 1 . Then, the above equation is written with grid-point indices for an arbitrary grid-point z p as:

    f k + p = f k + f Δ ¯ p , k + f Δ ¯ p , k 2 2 ! + f ( 3 ) Δ ¯ p , k 3 3 ! + f ( 4 ) Δ ¯ p , k 4 4 ! + R ,

    where the subscript p is an integer, Δ ¯ p , k = z p z k , and R means higher-order terms than the fourth.

    Second-order FDM in x-direction is approximated with a centered difference as following

    [ u h x ] i , k = u i , k h i + 1 , k h i 1 k 2 Ξ ,

    where (i, k) means the grid-index corresponding to (xi, zk) and Ξ implies the gridsize in the x-direction. Eq. (7) was obtained by neglecting the second-order or higher order terms from the Taylor expansions. Staggered-grid FDM on 3-points stencil for the vertical advection term is approximated in terms of an average of finite difference around the target point (Arakawa and Lamb, 1977; Haltiner and Williams, 1980; Tomita and Satoh, 2004; Cheong, 2006; Haiden et al., 2017):

    [ w h z ] i , k + 1 / 2 = 1 2 { [ w h z ] i , k + 1 + [ w h z ] i , k } = 0.5 w i , k + 1 ( h i , k + 1 h i , k Δ k + 1 ) + 0.5 w i , k ( h i , k h i , k 1 Δ k ) .

    As is expected from the Taylor series, this gives the second order accuracy for a uniform grid system and constant vertical velocity, while it is only first-orderaccuracy for nonuniform grid system regardless of the specific value of vertical velocity. Since the discretization error for the nonuniform grid is first-order, this scheme invites numerical diffusion in the vertical direction depending on the sign of advection speed. For uniform grid and constant speed (c), the finite difference form of the vertical advection term is written as

    c ( h i , k + 1 h i , k ) Δ k + 1 + Δ k + c ( h i , k h i , k 1 ) Δ k + Δ k 1 = c h + c 8 ( Δ k + 1 Δ k 1 ) h + R .

    When the above identity is applied to the advection equation, it is clear that the second-order derivative term acts as diffusion for c < 0 and upward increasing grid-size (i.e., the case of Δ k + 1 > Δ k 1 ) and c > 0 and downward increasing grid-size (i.e., the case of Δ k + 1 < Δ k 1 ). Negative diffusion is unavoidable in this scheme because three points including the target point at the center should be incorporated unlike the upwind differencing which uses only two points selectively.

    In the case where the advection velocity is not a constant value, Eq. (9) is rewritten as

    w i , k + 1 ( h i , k + 1 h i , k ) Δ k + 1 + Δ k + w i , k ( h i , k h i , k 1 ) Δ k + Δ k 1 = E 1 h + E 2 h + R ,


    E 1 = w i , k + 1 + w i , k 2
    E 2 = w i , k + 1 8 ( Δ k + 1 + Δ k ) w i , k 8 ( Δ k + Δ k 1 ) = Δ k 8 [ ( 1 + γ ) w i , k + 1 ( 1 + γ 1 ) w i , k ] = Δ k ( 1 + γ ) 8 [ w i , k + 1 1 γ w i , k ] .

    Equation (11) states that the coefficient of the discarded second-order term depends on the gridinterval and the advection velocity, which can be either positive (numerical adverse diffusion) or negative (numerical diffusion) as in the case of constant speed. This combined effect is expected to result in different behavior of the discretization error from the case of constant advection velocity.

    Initial Condition, Analytic Solution, and Model Parameters

    In numerical modeling of the atmospheric circulation, analytic solutions are often used to assess the accuracy of discretization (Browning et al., 1989; Williamson et al., 1992; Cheong, 2006; Cheong and Park, 2007; Flyer and Wright, 2007). In this study, the velocity field is specified as a pair of circulation cells of different sign (cell L and cell R), of which streamfunction, x-direction and z-direction velocity components are given as

    ψ ( x , z ) = { ψ L + ψ R ( 0 r 1 , r 2 0.5 H ) 0 ( r 1 , r 2 > 0.5 H ) }
    ψ L = + 2 + cos ( 2 π r 1 ) 2 ψ R = 2 + cos ( 2 π r 2 ) 2 r 1 = { ( x 0.5 H ) 2 + ( z 0.5 H ) 2 } r 2 = { ( x 1.5 H ) 2 + ( z 0.5 H ) 2 }


    u ( x , z ) ψ z = + sin ( 2 π r 1 ) π ( z 0.5 H ) r 1 sin ( 2 π r 2 ) π ( z 0.5 H ) r 2
    w ( x , z ) ψ x = + sin ( 2 π r 1 ) π ( x 0.5 H ) r 1 + sin ( 2 π r 2 ) π ( z 1.5 H ) r 2

    The flow pattern is very similar to the rotational flow on the sphere (Nair et al., 2005; Flyer and Wright, 2007; Cheong et al., 2015). The angular velocity of the advection flow is then written as

    ω ( r ) υ θ ( x , z ) r ψ r r = { π r 1 sin ( 2 π r 1 )    for cell L π r 2 sin ( 2 π r 2 )    for cell R

    Scalar field is given in the form of a Gaussian bell, which is expressed as

    h ( x , z , t ) = exp [ ( r r 0 ) 2 0.1 H ] × exp [ θ θ c ( r , t ) 2 D ( θ ) ] ,

    where r 0 = 3.0 H and D θ = ( 0.3 H ) 2 ( r 0.3 H ) 2 + 10 4 H , and θ c ( r , t ) denotes the azimuthal location of the Gaussian bell. Due to the circularity of the advection flow, it can be written as θ c ( r , t ) = ω ( r ) t , which enables to specify the scalar field at any time (ω=angular speed).

    The horizontal domain is set to be twice of the vertical domain so that two circular advection flows can be included in the model domain. The model equation is scaled is such a way that both the maximum advection flow and the vertical range are unity. In the control run, the number of the grid-point in z-direction is set Nz= 30 and that of x-direction is set as Nx= 2Nz, with which the spatial resolutions in both directions are the same for the uniform grid spacing. Time integration is carried out until the maximum advection flow transports the scalar field along azimuthal direction by the nondimensional distance of 0.3 at the radial distance of r = 0.25H. The time-step size is given large enough to make the CFL number be unity for the smallest grid-size, which is smaller than the marginal value ( 3 ) for the stability (Nair et al., 2005; Kang and Cheong, 2017).

    Results and Discussion

    As depicted in section 3, the numerical error will depend on grid structure (that is, non-uniformity) as well as the advection flow configuration even though the discretization method is the same. This issue was taken into consideration by defining the initial scalar field at various locations. Figure 3 illustrates two examples of the initial distribution of the scalar field. In the first example, which corresponds to the control run, two scalar fields are located in the same position and the scalar variable is advected away from the boundary, as indicated by arrows. In the case of the nonuniform grid, this means that it is transported from fine vertical-resolution region to coarse verticalresolution region. In the lower-panel example, the left cell has the initial scalar field in the same position the upper panel. On the other hand, the right cell has the scalar field near the upper boundary so that it is transported from the coarse vertical-resolution region to the fine vertical-resolution region. It is expected that transport of the scalar field across different resolution grids will result in different error behavior. To assess the accuracy, the root mean squared error (RMSE) is measured for each circulation cell:

    RMSE = i = 1 N x k = 1 N z ( h s ( i , k ) h r ( i , k ) N x N z ) 2 ,

    where hs and hr mean the simulated field and reference field, respectively.

    Figure 4 presents the simulation results for the control run. As is expected from flow configuration, the error distribution is exactly symmetric with respect to the vertical line x = L/2 (=0.5 in non-dimensional unit). This is because the centered FDM is employed in x-direction, which, unlike the vertical discretization, has the same error for both directions (negative and positive x-direction). Summary of simulations with various initial location of tracer, spatial resolutions, and non-uniformity, is shown in Table 1. The results clearly exhibit the second-order error convergence for the uniform grid. For instance, the error of case III is reduced by a factor of about 4 in comparison with the case II. The same holds true for the cases of III and IV. It is noticeable that the error for the nonuniform grid is quite sensitive to the initial location of the tracer as can be seen the results of the cases from V to X. In the cases of the nonuniform grid, the error of downdraft motion in the stretching-grid region (i.e., the grid-size increases with height) is larger than the downdraft in the shrinking-grid region (i.e., the gridsize decreases with height). This tendency becomes more conspicuous for increased nonuniformity, e.g., the cases of from V to X. Such a behavior is well illustrated in Fig. 5 which corresponds to the case V. Error patterns appear in wave train with alternating sign, and more contour lines, being indicative of larger error, is clearly seen in the left circulation cell. The reason for this asymmetric error behavior may be that the initial condition of Gaussian bell is resolved not in an identical accuracy. It is certain that the initial tracer located in the fine-grid region better resolves the Gaussian bell than that located in the coarse-grid region. This is well illustrated in Fig. 6 where the results of case V are plotted. Obviously, the initial tracer over the left cell is defined on smaller number of grid-points than the right cell.

    To identify the combined effect of nonuniformity of grid and the stretching or shrinking motion in the vertical direction more clearly, an experiment with circular (i.e., axisymmetric) pattern of tracer distribution was conducted. Since both the flow pattern and the tracer distribution are axisymmetric, the error pattern also should be axisymmetric if the error source associated with the spatial discretization is axisym- metric. Result of this case is shown in Fig. 7, where the error pattern is not found to be axisymmetric. This is attributed to the non-uniformity of the grid and stretching or shrinking motion in the vertical direction. Combined effect of these two factors is represented with Eq. (11b). The error source may be termed as either numerical diffusion (E2< 0) or numerical adverse-diffusion (E2> 0). The largest error, therefore, is directly associated with the largest value of E2, which may be observed in the updraft or downdraft near the boundaries. The adverse diffusion causes the time integration unstable while the diffusion does not. Of the two factors in Eq. (11b), the stretching or shrinking of the flow in the vertical direction is more important than the non-uniformity of the grid as for the two dimensional circulation. This is because the ratio of stretching velocity, for adjacent two grid points, is well above 3 or 4 for in the level near the boundaries, being much larger than the value of γ. Since the error pattern is advected along the circular trajectory, it can be imagined easily that the error at the initial stage is distributed near the boundaries. That is, the vertical distribution of |E2| is severely confined to the region of strong velocity-stretching. If the results are applied to the atmospheric circulation model, it may be stated that the updraft in the boundary layer accompanied by the horizontal convergence is vulnerable to large numerical error.

    Concluding Remarks

    The effect of nonuniform grid on the two dimensional transport phenomenon was investigated based on the numerical simulations and error analysis. Main concern of the study was the vertical discretization and nonuniform staggered-grid which resembled the weather forecast model in a qualitative sense. Discretization method in time was the thirdorder Runge-Kutta scheme, while for spatial discretization the FDM on the staggered-grid 3-points stencil in the vertical direction and the centered FDM in the horizontal direction were used. Error analysis indicated that the vertical discretization produces numerical adverse-diffusion or diffusion when the vertical-flow slope (that is, vertical convergence or divergence) is large. This effect tends to be strengthened by the presence of non-uniformity of the grid. However, the effect of the latter was smaller than that of the effect of vertical-flow slope for the typical non-uniformity of the grid. This result implies that in the numerical model the rising or sinking motion, accompanied by horizontal convergence or divergence, respectively, near the boundary, may be vulnerable to the numerical instability due to the adverse diffusion. And the possibility of the instability becomes large in the case of non-uniformity of grid-size.


    This work was supported by a Research Grant of Pukyong National University (year 2017-2018). Anonymous reviewers are acknowledged for their useful comments.



    Nonuniform grid system of the present study: The grid interval in z-direction is nonuniform with decreasing toward the upper and lower boundaries while the interval in x-direction is uniform.


    Location of variables. Cross mark denotes the scalar variable, and open circle and solid circle denote x-direction velocity and z-direction velocity, respectively.


    Examples of initial tracer distribution (shaded circle), the streamfunction (solid and broken lines), and corresponding advection direction (thick arrow). See the text for more detailed explanation.


    Simulation results of the reference case. Contour plot superimposed on the velocity vector field represents the initial distribution of the tracer. Contour plot in the right bottom panel is the error of the tracer field, where positive (negative) values are in solid (dashed) lines.


    Same as Fig. 4 except for the case II (Nx= 60, Nz= 30, γ = 1.00).


    Same as Fig. 4 except for the case V (Nx= 60, Nz= 30, γ = 1.05).


    Same as the case VII (Nx= 120, Nz= 60, γ = 1.1) except for the circular pattern of the initial condition.


    Cases of numerical experiments where Nx (Nz) is number of horizontal (vertical) grid and γ means the nonuniformity of vertical grid (for definition, see the text). The CFL number is given as unity for all experiments. Symbol H (L) implies the size of horizontal (vertical) domain


    1. Arakawa, A. , 1972, Design of the UCLA general circulation model. UCLA Meteorology Department Tech. Rep. 7, 116 pp.
    2. Arakawa, A. , and Lamb, V.R. , 1977: Computational design of the basic dynamical processes of the UCLA general circulation model. Methods in Computational Physics,Vol. 17, J. Chang, Ed., Academic Press, 173-265.
    3. Browning, G.L. , Hack, J.J. , and Swarztrauber, P.N. , 1989,  A comparison of three numerical methods for solving differential equations on the sphere . Monthly WeatherReview, 117, 1058-1075.
    4. Cheong, H.B. , 2006: A dynamical core with double Fourier series: Comparison with the spherical harmonics method . Monthly Weather Review, 134, 1299-1315.
    5. Cheong, H.B. , Kong, H.J. , Kang, H.G. , and Lee, J.D. , 2015, Fourier finite-element method with linear basis functions on a sphere: Application to elliptic and transport equations . Monthly Weather Review, 143,1275-1294.
    6. Cheong, H.B. and Park, J.R. , 2007, Geopotential field in nonlinear balance with the sectoral mode of Rossby- Haurwitz wave on the inclined rotation axis . Journal of the Korean Earth Science Society, 28, 936-946.
    7. Choi, S.J , and Hong, S.Y. , 2016, A global non-hydrostatic dynamical core using the spectral element method on a cubed-sphere grid . Asia Pacific Journal of AtmosphericSciences, 52, 291-307.
    8. Durran, D.R. , 1999, Numerical methods for wave equations in geophysical fluid dynamics. Springer, 465pp.
    9. Flyer, N. , and Wright, G.B. , 2007, Transport schemes on a sphere using radial basis functions . J. Computational Physics, 226, 1059-1084.
    10. Gal-Chen, T. , and Somerville, R. , 1975: On the use of a coordinate transformation for the solution of the Navier- Stokes equations . Journal of Computational Physics, 17, 209-228.
    11. Gates, W.L. , 1992, AMIP: The Atmospheric Model Intercomparison Project . Bulletin of American Meteorological Society, 73, 1962-1970.
    12. Haiden, T. , Janousek, M. , Bidlot, J.R. , Ferranti, L. , Prates, F. , Vitart, F. , Bauer, P. , and Richardson, D. , 2017, Evaluation of ECMWF forecasts, including 2016-2017 upgrades. Technical Memorandum 817, ECMWF, 58 pp.
    13. Haltiner, G.J. , and Williams, R.T. , 1980, Numerical weather prediction and dynamic meteorology. John Wiley and Sons, 477 pp.
    14. Hoskins, B.J. , and Simmons, A.J. , 1975, A multilayer spectral model and the semiimplicit method . Quarterly Journal of Royal Meteorological Society, 101, 637-655.
    15. Kang, H.G. , and Cheong, H.B. , 2017, An efficient implementation of a high-order filter for a cubed-sphere spectral element model . J. Computational Physics, 332, 66-82.
    16. Kang, H.G. , and Cheong, H.B. , 2018, Effect of high-order filter on a cubed-sphere spectral element dynamicalcore . Monthly Weather Review, 146, 2047-2064.
    17. Klemp, J.B. , 2011, A Terrain-Following Coordinate with Smoothed Coordinate Surfaces . Monthly WeatherReview, 139, 2163-2169.
    18. Leuenberger, D. , Koller, M. , Fuhrer, O. , and Schr, C. , 2010, A generalization of the SLEVE vertical coordinate . Monthly Weather Review, 138, 3683-3689.
    19. Marras, S.M. , Kopera, A. , and Giraldo, F.X. , 2015, Simulation of shallow-water jets with a unified elementbased continuous/discontinuous Galerkin model with grid flexibility on the sphere . Quarterly Journal of Royal Meteorological Society, 141, 1727-1739.
    20. Nair, R.D. , Thomas, S.J. , and Loft, R.D. , 2005, A discontinuous Galerkin transport scheme on the cubedsphere . Monthly Weather Review, 133, 814-828.
    21. Rivier, L. , Loft, R. , and Polvani, L.M. , 2002, An efficient spectral dynamical core for distributed memory computers . Monthly Weather Review, 130, 1384-1396.
    22. Saito, K. , Fujita, T. , Yamada, Y. , Ishida, J.I. , Kumagai, Y. , Aranami, K. , Ohmori, S. , Nagasawa, R. , Kumagai, S. , Muroi, C. , Kato, T. , Eito, H. , and Yamazaki, Y. , 2006, The operational JMA nonhydrostatic mesoscale model . Monthly Weather Review, 134, 1266-1298.
    23. Simmons, A.J. , and Burridge, D.M. , 1981, An energy and angular-momentum conserving vertical finite-difference scheme and hybrid vertical coordinates . Monthly Weather Review, 109, 758-766.
    24. Skamarock, W.C. , and Klemp, J.B. , 2008, A time-split nonhydrostatic atmospheric model for weather research and forecasting applications . J. Computational Physics, 227, 3465-3485.
    25. Smagorinsky, J. , 1963, General circulation experiments with the primitive equations: I. The basic experiment . Monthly Weather Review, 91, 99-164.
    26. Tomita, H. , and Satoh, M. , 2004, A new dynamical framework of nonhydrostatic global model using the icosahedral grid . Fluid Dynamics Research, 34, 357-400.
    27. Williamson, D.L. , Drake, J.B. , Hack, J.J. , Jakob, R. , and Swarztrauber, P.N. , 1992, A standard test set fornumerical approximations to the shallow water equations in spherical geometry . J. ComputationalPhysics, 102, 211-224.
    28. Zhang, H. , and Rancic, M. , 2007, A global eta model on quasi-uniform grids . Quarterly Journal of RoyalMeteorological Society, 133, 517-528.