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

LI XiaohanPENG Xindongand LI Xingliang

1State Key Laboratory of Severe Weather,Chinese Academy of Meteorological Sciences,Beijing100081

2Center of Numerical Weather Prediction,China Meteorological Administration,Beijing100081

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

LI Xiaohan1,PENG Xindong?1,and LI Xingliang1,2

1State Key Laboratory of Severe Weather,Chinese Academy of Meteorological Sciences,Beijing100081

2Center of Numerical Weather Prediction,China Meteorological Administration,Beijing100081

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 ef fi ciency.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 coef fi cients of the equation are quite different from those in the original model,the computational procedureof the newcore is just thesame.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 fl ow and development of both the mountain-induced Rossby wave and Rossby–Haurwitz wave con fi rms the appropriate installation of the 3D Coriolis terms in the semi-implicit semi-Lagrangian dynamic core on the Yin-Yang grid.

Yin-Yang grid,semi-implicit semi-Lagrangian,nonhydrostatic,dynamic core

1. Introduction

Thedynamiccoreistheheartoftheatmosphericmodel.It determines the computing characteristics such as the numerical accuracy and computational ef fi ciency.Use of a highorder spatial fi nite differencing scheme,semi-Lagrangian transport,advanced temporal integration,and other stateof-the-art techniques has improved model representation signi fi cantly.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 nonhydrostatic 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 boundaryarrangement.TheGRAPESdynamiccoreissolved with the semi-implicit semi-Lagrangian(SISL)method containing the spherical departure-point determination(Ritchie and Beaudoin,1994)and the non-centered fi nite 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 dif fi culties,such as pole singularity and the conver-gence 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 whendealingwiththesphericaldeparture-pointintheRitchie 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.Themostfamiliardesignsofquasi-uniformgrids 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 re fi nement 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 computingprocedureoftheSISL methodwhen a3DCoriolis forcing is added.Rede fi nition of the coef fi cients 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 coef fi cients,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, fi rst 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 fi t 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 fi nite 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 de fi nitely 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- fl ux balance between the component zones.Peng et al.(2006)developed a numerical constraint to guarantee that the fl uxes 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 signi fi cantly 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 de fi ned for the cubic Lagrange interpolation.The coordinate conversion between Yin and Yang grids can be expressed as

whereλand?represent longitude and latitude,respectively, and the subscriptseandodenote 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

and vice versa.

3. Dynamical frame of the GRAPESYY

3.1.Prognostic equations

The nonhydrostatic governing equations of the atmosphere on a sphere read

on a spherical coordinate system,where Π is the Exner function,θis the potential temperature,cp=1004.64 J kg-1K-1represents the speci fi c heat at constant pressure,u,vare horizontal winds andwis the vertical speed,tmeans time,g=9.80616 m s-2is the gravitational acceleration,ris the radius vector of the spherical coordinate,D3denotes the 3D divergence,fis Coriolis parameter and

Rdrepresents the ideal gas constant for dry air,QTshows the source or sink term of heat andFTis 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 fl exibility of the dynamic core.After discretization with the SISL method,the dynamical equations for the prognostic variablesu,v,w,Π′andθ′at time leveln+1 are written as

in a height-based terrain-following coordinatewhere Π′(θ′)is a perturbation of the Exner function(potential temperature)from its reference state?Π(?θ).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 heightdependent function of each variable,which is the horizontal average of the given initial fi elds and is different from that in the original model.In the equations,?tis the time step andαεthe contribution adjustment factor.Information on the departure point at time levelnand the nonlinear terms are included inAx(x=u,v,w,Π′,θ′),which remains the same as intheoriginalGRAPESmodel(XueandChen,2008).Linear termsLx(x=u,v,w,Π′,θ′)at time leveln+1 are

where

is a residual fraction of the reference atmosphere from the hydrostatic balance state,and

φsxandφsyare the topographic slope,zTis the model top,andzsissurfaceheight.Thetermscontainingf?andfλarenewly introduced into the equations in comparison to the GRAPES model.Considering Eqs.(13)–(15),the linear equation set of Eqs.(8)–(10)can be solved accordingly:

where

and

δ=αε?tis de fi ned,and matrixis given as

and

wherezis the height level.Substitutewinto the thermodynamic equation[Eq.(17)]to obtainθ′at the next time level:

where

Whenu,v,wandLΠin Eq.(11)are substituted with Eqs. (18)–(20)and(16),a Helmholtz equation is deduced as

where

The related 19 Π grid points,which serves the numerical solution of Eq.(39),are displayed in Fig.1.It is worth noting that Eqs.(21)–(43)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 incompleteLL Ufactorization(IL U)(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-15for the numerical convergence measurement.

3.2.Semi-Lagrangian advection

3.2.1.Semi-Lagrangian transport on the Yin-Yang grid

Both the nonlinear terms and the departure-point-related termsareincludedinAxinEqs.(8)–(12).Thedeparturepoint is calculated according to Ritchie and Beaudoin(1994).Halo regions are de fi ned for each component zone to avoid multitime 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,rm)at the half-time level on the sphere considering

where(λa,?a,ra)represents the arrival point.

(2)Determine the velocity components at the midpoint (um,vm,wm)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,rd)is then de fi ned as

Note that steps(1)and(3)are the same as in the original model,while step(2)is modi fi ed 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 de fi ne the coef fi cients 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 calculateAu,v,won the Yin-Yang grid.For scalars,the termAθ,Πcan be interpolated at the departure point(denoted with superscript“*”)directly:

whereLx(x=θ,Π)is the linear term,Nx(x=θ,Π)is the nonlinear variation and

We note that the termsξu0,v0,w0all containAu,v,w,θin Eqs. (24),(28)and(32).The prediction of 3D momentum calls for a good description ofAu,v,w,θ,which depends on the computation of

at their departure points，denote by.

In the non-uniform vertical direction,a cubic Lagrangian interpolation,

is used to ensure high-order accuracy.A second-order formulation for the vertical gradient of Π is given as

This ensures second-order accuracy without guaranteeing quantity conservation in the vertical transport.An alternative choice that achieves exact conservation is

but the accuracy decreases for non-uniform grids,where the Πkis located at the middle level of?zkand?zk+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 modi fi cation 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 de fi ned 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.

4.1.Steady-state geostrophic fl ow

This test is a 3D extension of Test 2 in Williamson et al. (1992).The initial state is de fi ned as

whereu0is set to be 20 m s-1andαis fl ow orientation angle, which is 0 here.Thirty-day integration results of Π′,u,v,wand the differences between the numerical solution and the exactone,withatimestepof1800s,areshowninFigs.3a–d. Corresponding results from GRAPES are presented in Figs. 3e–h for comparison.The Π′in Fig.3a shows its zonal parallel contours after the 30-day integration,while the absolute error is about-1.0×10-4at the equator and 5.0×10-5in midlatitude regions.Zonal wind(Fig.3b)keeps its initial state with a perturbation of-0.12 m s-1at 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-3and 10-7m s-1,respectively.Meanwhile,errors of meridional and vertical winds reach 0.1 and 2×10-3m s-1in 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 signi fi cant.Consequently, a boundary trail is revealed in Figs.3c and d.The time series of the corresponding error norms?1,?2and?∞of both the scalar and vectors are given for GRAPES YY in Fig.4.Error norms increase with time.The?1and?2norms of the Π′are 0.002 and 0.0025 at day 30,while those of the velocity are 0.009 and 0.0095,respectively.This numerical test con fi rms the stability of the SISL method and the proper installation of the 3D Coriolis force in the nonhydrostatic frame.

In the zonal fl ow case,small interpolation error at the boundaries is clearly displayed with the meridional wind and vertical velocity.We fi nd the error ofvandwto be 10-3and 10-7,respectively,which is negligible in comparison with theucomponent.But does the error destroy the model stability in a non-zero meridional wind case?We also show the numerical results of the balanced fl ow with an angle of 45°in Fig.5.This con fi guration 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 Π′,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 fi elds.

4.2.Zonal fl ow over a mountain

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

whereh0=2000 m determines the peak height of the mountain andR=1500 km is the mountain half width;D,thedistance to the mountain center(λc,?c)=(π/2,π/6):

The pressure fi eld is initially given a hydrostatic balance state,

whereN=0.0182 s-1is the Brunt–Valsala frequency,κ= 2/7,u0=20 m s-1,a=6371.229 km is mean radius of earth,? is earth’s angular velocity andpsp=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 componentsuandvat 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 lowresolution con fi guration is used here.The fi gures 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 con fi guration.Therefore,we can again con fi rm good numerical performance via this test case.

4.3.Three-dimensional Rossby–Haurwitz wave

The initial velocity fi eld is given by

wherec=4 denotes the wave number,andM=K≈1.962×10-6s-1.Pro fi les of the temperature and pressure are given by

wherepref=955 hPa andT0=288 K with the moistadiabatic lapse rate Γ=0.0065 km-1.The geopotential is given as

where

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,andvat 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 fi eld 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.

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 ef fi ciency.The iteration before its convergence is listed in Fig.8 for the fi rst 100-step integration.Rapid convergence of the solvers is achieved for the overset grid,and the iteration tends to decrease with time.The test of the zonal fl ow over a mountain shows more iterations than the others for its strong time-dependent current.

5. 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 fi xed 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 ef fi cient 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 coef fi cients fortheHelmholtzequationintheclassicalSchwarzmethodis 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.

Acknowledgements.The authors appreciate the two reviewers’constructive comments and valuable suggestions.This study was supported by the National Natural Science Foundation of China (Grant No.41175095),the National Key Technology R&D Program(Grant No.2012BAC22B01)and a research project of the Chinese Academy of Meteorological Sciences(Grant No.2014Z001).

REFERENCES

Arakawa,A.,and 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.,and N.A.Philips,1953:Numerical integration of the quasi-geostrophic equations for barotropic and simple baroclinic Flow.J.Meteor.,10,71–99.

Chen,D.,and Coauthors,2008:New generation of multi-scale NWP system(GRAPES):General scienti fi c design.Chinese Sci.Bull.,53(22),3433–3445.

Jablonowski,C.,P.Lauritzen,R.Nair,and M.Taylor,2008:Idealized test cases for the dynamical cores of atmospheric general circulationmodels:AproposalfortheNCARASP2008summer colloquium.[Available online at http://esse.engin.umich. edu/admg/publications.php.]

Kageyama,A.,and 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 fi nite-volume shallow-water model on the Yin-Yang overset spherical grid.Mon.Wea.Rev.,136, 3066–3086.

Liu,Y.,and J.W.Cao,2008:ILU preconditioner for NWP system:GRAPES.Computer Engineering and Design,29(3), 731–734.(in Chinese)

McDonald,A.,and 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?ot′e,and M.Gander,2008:Optimized Schwarz methods with an overset grid for shallowwater equations:Preliminary results.Appl.Numer.Math., 58(4),459–471.

Ritchie,H.,and 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 fi nite-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.,and 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.,and D.H.Chen,2008:Design of GRAPES dynamical frame and the experiments.Scienti fi c Design and Application of Numerical Prediction System GRAPES,J.H.Wang,Ed., Science Press,Beijing,65–136.(in Chinese)

:Li,X.H.,X.D.Peng,and X.L.Li,2015:An improved dynamic core for GRAPES on the Yin-Yang grid.Adv. Atmos.Sci.,32(5),648–658,

10.1007/s00376-014-4120-5.

(Received 11 June 2014;revised 22 September 2014;accepted 26 September 2014)

?Corresponding author:PENG Xindong

Email:pengxd@cams.cma.gov.cn