Advanced Search
Article Contents

Improvement of the Semi-Lagrangian Advection Scheme in the GRAPES Model: Theoretical Analysis and Idealized Tests

Fund Project:

doi: 10.1007/s00376-013-3086-z

  • The Global/Regional Assimilation and PrEdiction System (GRAPES) is the new-generation numerical weather prediction (NWP) system developed by the China Meteorological Administration. It is a fully compressible non-hydrostatical global/regional unified model that uses a traditional semi-Lagrangian advection scheme with cubic Lagrangian interpolation (referred to as the SL_CL scheme). The SL_CL scheme has been used in many operational NWP models, but there are still some deficiencies, such as the damping effects due to the interpolation and the relatively low accuracy. Based on \hboxReich's semi-Lagrangian advection scheme (referred to as the R2007 scheme), the Re_R2007 scheme that uses the low- and high-order B-spline function for interpolation at the departure point, is developed in this paper. One- and two-dimensional idealized tests in the rectangular coordinate system with uniform grid cells were conducted to compare the Re_R2007 scheme and the SL_CL scheme. The numerical results showed that: (2) the damping effects were remarkably reduced with the Re_R2007 scheme; and (3) the normalized errors of the Re_R2007 scheme were about 7.5 and 3 times smaller than those of the SL_CL scheme in one- and two-dimensional tests, respectively, indicating the higher accuracy of the Re_R2007 scheme. Furthermore, two solid-body rotation tests were conducted in the latitude-longitude spherical coordinate system with non-uniform grid cells, which also verified the Re_R2007 scheme's advantages. Finally, in comparison with other global advection schemes, the Re_R2007 scheme was competitive in terms of accuracy and flow independence. An encouraging possibility for the application of the Re_R2007 scheme to the GRAPES model is provided.
    摘要: The Global/Regional Assimilation and PrEdiction System (GRAPES) is the new-generation numerical weather prediction (NWP) system developed by the China Meteorological Administration. It is a fully compressible non-hydrostatical global/regional unified model that uses a traditional semi-Lagrangian advection scheme with cubic Lagrangian interpolation (referred to as the SL_CL scheme). The SL_CL scheme has been used in many operational NWP models, but there are still some deficiencies, such as the damping effects due to the interpolation and the relatively low accuracy. Based on \hboxReich's semi-Lagrangian advection scheme (referred to as the R2007 scheme), the Re_R2007 scheme that uses the low- and high-order B-spline function for interpolation at the departure point, is developed in this paper. One- and two-dimensional idealized tests in the rectangular coordinate system with uniform grid cells were conducted to compare the Re_R2007 scheme and the SL_CL scheme. The numerical results showed that: (2) the damping effects were remarkably reduced with the Re_R2007 scheme; and (3) the normalized errors of the Re_R2007 scheme were about 7.5 and 3 times smaller than those of the SL_CL scheme in one- and two-dimensional tests, respectively, indicating the higher accuracy of the Re_R2007 scheme. Furthermore, two solid-body rotation tests were conducted in the latitude-longitude spherical coordinate system with non-uniform grid cells, which also verified the Re_R2007 scheme's advantages. Finally, in comparison with other global advection schemes, the Re_R2007 scheme was competitive in terms of accuracy and flow independence. An encouraging possibility for the application of the Re_R2007 scheme to the GRAPES model is provided.
  • 加载中
  • Bermejo, R., and A. Staniforth, 1992: The conversion of semi-Lagrangian advection schemes to quasi-monotone schemes. Mon. Wea. Rev., 120(12),2622-2632.
    Chen, D. H., and Coauthors, 2008: New generation of multi-scale NWP system (GRAPES): General scientific design. Chinese Sci. Bull., 53(23),3433-3445.
    Cotter, C., J. Frank and S. Reich, 2007: The remapped particlemesh semiLagrangian advection scheme. Quart. J. Roy. Meteor. Soc., 133(622),251-260.
    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(6),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 shallowwater model on the Yin-Yang overset spherical grid. Mon. Wea. Rev., 136(9),3066-3086.
    Lin, S. J., and R. B. Rood, 1996: Multidimensional flux-form semi-Lagrangian transport schemes. Mon. Wea. Rev., 124(10),2046-2070.
    Liu, H. Y., J. S. Xue, J. F. Gu, and H. M. Xu, 2012: Radar data assimilation of the GRAPES model and experimental results in a typhoon case. Adv. Atmos. Sci., 29(3),344-358, doi: 10.1007/s00376-011-1063-y.
    McDonald, A., 1984: Accuracy of multiply-upstream, semi-Lagrangian advective schemes. Mon. Wea. Rev, 112, 1267-1275.
    McDonald, A., 1987: Accuracy of multiply-upstream semi-Lagrangian advective schemes II. Mon. Wea. Rev, 115, 1446-1275.
    McDonald, A., and J. R. Bates, 1989: Semi-Lagrangian integration of a gridpoint shallow water model on the sphere. Mon. Wea. Rev, 117(2),130-137.
    Nair, R. D., and B. Machenhauer, 2002: The mass conservative cell-integrated semi-Lagrangian advection scheme on the sphere. Mon. Wea. Rev., 130(4),649-667.
    Powell, M. J. D., 1981: Approximation Theory and Methods. Cambridge University Press, 339 pp.
    Rasch, P. J., 1994: Conservative shape-preserving two-dimensional transport on a spherical reduced grid. Mon. Wea. Rev., 122(7),1337-1550.
    Reich, S., 2007: An explicit and conservative remapping strategy for semi-Lagrangian advection. Atmospheric Science Letters, 8(3),58-63.
    Ritchie, H., and C. Beaudoin, 1994: Approximations and sensitivity experiments with a baroclinic semi-Lagrangian spectral model. Mon. Wea. Rev., 122(11),2391-2399.
    Robert, A., 1981: A stable numerical integration scheme for the primitive meteorological equations. Atmos.-Ocean, 19, 35-46.
    Robert, A., 1982: Semi-Lagrangian and semi-implicit numerical integration scheme for the primitive meteorological equation. Quart. J. Roy. Meteor. Soc., 60, 319-324.
    Staniforth, A., and J. C#cod#x0020f;t#cod#x000e9;, 1991: Semi-Lagrangian integration schemes for atmospheric models——A review. Mon. Wea. Rev., 119, 2206-2223.
    Williamson, D. L., J. B. Drake, J. J. Hack, R. Jackob, and P. Swarztrauber, 1992: A standard test set for numerical approximations to the shallow water equations in spherical geometry. J. Comput. Phys., 102(2),211-224.
    Xue, J. S., 2004: Progresses of researches on numerical weather prediction in China: 1999-2002. Adv. Atmos. Sci., 21(4),467-474.
    Xue, J. S., and Y. Liu, 2007: Numerical weather prediction in China in the new century progress, problems and prospects. Adv. Atmos. Sci., 24(7),1099-1108, doi: 10.1007/s00376-007-1099-1.
    Yang, L. J., and X. S. Shen, 2011: The construction of SCM in GRAPES and its applications in two field experiment simulations. Adv. Atmos. Sci., 28(4),534-550, doi: 10.1007/s00376-010-0062-8.
    Zerroukat, M., N. Wood and A. Staniforth, 2004: SLICE-S: A semi-Lagrangian inherently conserving and efficient scheme for transport problems on the sphere. Quart. J. Roy. Meteor. Soc., 130(602),2649-2664.
  • [1] Yi ZHANG, Rucong YU, Jian LI, 2017: Implementation of a Conservative Two-step Shape-Preserving Advection Scheme on a Spherical Icosahedral Hexagonal Geodesic Grid, ADVANCES IN ATMOSPHERIC SCIENCES, 34, 411-427.  doi: 10.1007/s00376-016-6097-8
    [2] ZHANG Kai, WAN Hui, WANG Bin, ZHANG Meigen, 2008: Consistency Problem with Tracer Advection in the Atmospheric Model GAMIL, ADVANCES IN ATMOSPHERIC SCIENCES, 25, 306-318.  doi: 10.1007/s00376-008-0306-z
    [3] Ji Zhongzhen, Wang Bin, 1997: Multispectrum Method and the Computation of Vapor Equation, ADVANCES IN ATMOSPHERIC SCIENCES, 14, 563-568.  doi: 10.1007/s00376-997-0074-1
    [4] Eric P. CHASSIGNET, Xiaobiao XU, 2021: On the Importance of High-Resolution in Large-Scale Ocean Models, ADVANCES IN ATMOSPHERIC SCIENCES, 38, 1621-1634.  doi: 10.1007/s00376-021-0385-7
    [5] 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
    [6] Fangyuan CHENG, Qinghua YANG, Changwei LIU, Bo HAN, Shijie PENG, Guanghua HAO, 2023: Evaluating Parameterizations for Turbulent Fluxes over the Landfast Sea-Ice Surface in Prydz Bay, Antarctica, ADVANCES IN ATMOSPHERIC SCIENCES, 40, 1816-1832.  doi: 10.1007/s00376-023-2299-z
    [7] 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
    [8] 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
    [9] Qian Yongfu, 1988: THE LOCAL SPLINE VERTICAL INTERPOLATION METHOD OF TEMPERATURE AND GEOPOTENTIAL HEIGHT FIELDS AND THE TIME-DEPENDENT DIFFERENCE FORM OF THE HYDROSTATIC EQUATION, ADVANCES IN ATMOSPHERIC SCIENCES, 5, 159-170.  doi: 10.1007/BF02656778
    [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] 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
    [12] LI Qian, SUN Shufen, DAI Qiudan, 2009: The Numerical Scheme Development of a Simplified Frozen Soil Model, ADVANCES IN ATMOSPHERIC SCIENCES, 26, 940-950.  doi: 10.1007/s00376-009-7174-z
    [13] S.K. Sinha, D.R. Talwalkar, S.G. Narkhedkar, S. Rajamani, 1989: A Scheme for Objective Analysis of Wind Field Incorporating Multi-Weighting Functions in the Optimum Interpolation Method, ADVANCES IN ATMOSPHERIC SCIENCES, 6, 435-446.  doi: 10.1007/BF03342547
    [14] Qian Yongfu, Zhong Zhong, 1986: GENERAL FORMS OF DYNAMIC EQUATIONS FOR ATMOSPHERE IN NUMERICAL MODELS WITH TOPOGRAPHY, ADVANCES IN ATMOSPHERIC SCIENCES, 3, 10-22.  doi: 10.1007/BF02680042
    [15] Xu Hong, Li Hongji, Wang Ronghua, 1989: A Numerical Method of Statistical Pattern Recognition, ADVANCES IN ATMOSPHERIC SCIENCES, 6, 483-492.  doi: 10.1007/BF02659082
    [16] 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
    [17] Huang Ronghui, Li Xu, Yuan Chongguang, Lu Riyu, Moon Sung-Euii, Kim Ung-Jun, 1998: Seasonal Prediction Experiments of the Summer Droughts and Floods during the Early 1990’s in East Asia with Numerical Models, ADVANCES IN ATMOSPHERIC SCIENCES, 15, 433-446.  doi: 10.1007/s00376-998-0025-5
    [18] Yu Rucong, 1994: A Two-Step Shape-Preserving Advection Scheme, ADVANCES IN ATMOSPHERIC SCIENCES, 11, 479-490.  doi: 10.1007/BF02658169
    [19] Zhao Li, Zhao Sixiong, 1995: Numerical Experiments of Meiyu(Baiu) Rainfall by Quasi-Lagrangian Limited Area Model with Terrain, ADVANCES IN ATMOSPHERIC SCIENCES, 12, 57-66.  doi: 10.1007/BF02661287
    [20] Yi Zengxin, T. Warn, 1987: A NUMERICAL METHOD FOR SOLVING THE EVOLUTION EQUATION OF SOLITARY ROSSBY WAVES ON A WEAK SHEAR, ADVANCES IN ATMOSPHERIC SCIENCES, 4, 43-54.  doi: 10.1007/BF02656660

Get Citation+

Export:  

Share Article

Manuscript History

Manuscript received: 01 May 2013
Manuscript revised: 10 July 2013
通讯作者: 陈斌, bchen63@163.com
  • 1. 

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

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

Improvement of the Semi-Lagrangian Advection Scheme in the GRAPES Model: Theoretical Analysis and Idealized Tests

    Corresponding author: CHEN Dehui, chendh@cma.gov.cn
  • 1. Chinese Academy of Meteorological Sciences, Beijing 100081 ; 
  • 2. Numerical Weather Prediction Center, China Meteorological Administration, Beijing 100081 ; 
  • 3. University of Chinese Academy of Sciences, Beijing 100049
Fund Project:  We thank the anonymous reviewers for their constructive comments and also gratefully acknowledge Prof. Sebastian REICH (currently in the Department of Mathematics at Potsdam University) for his kind support and valuable suggestions. This work was jointly sponsored by the Key Project of the Chinese National Programs for Fundamental Research and Development (973 Program; Grant No. 2013CB430106), the Key Project of the Chinese National Science Technology Pillar Program during the Twelfth Five-year Plan Period (Grant No. 2012BAC22B01), and the National Natural Science Foundation of China ( Grant No. 41375108).

Abstract: The Global/Regional Assimilation and PrEdiction System (GRAPES) is the new-generation numerical weather prediction (NWP) system developed by the China Meteorological Administration. It is a fully compressible non-hydrostatical global/regional unified model that uses a traditional semi-Lagrangian advection scheme with cubic Lagrangian interpolation (referred to as the SL_CL scheme). The SL_CL scheme has been used in many operational NWP models, but there are still some deficiencies, such as the damping effects due to the interpolation and the relatively low accuracy. Based on \hboxReich's semi-Lagrangian advection scheme (referred to as the R2007 scheme), the Re_R2007 scheme that uses the low- and high-order B-spline function for interpolation at the departure point, is developed in this paper. One- and two-dimensional idealized tests in the rectangular coordinate system with uniform grid cells were conducted to compare the Re_R2007 scheme and the SL_CL scheme. The numerical results showed that: (2) the damping effects were remarkably reduced with the Re_R2007 scheme; and (3) the normalized errors of the Re_R2007 scheme were about 7.5 and 3 times smaller than those of the SL_CL scheme in one- and two-dimensional tests, respectively, indicating the higher accuracy of the Re_R2007 scheme. Furthermore, two solid-body rotation tests were conducted in the latitude-longitude spherical coordinate system with non-uniform grid cells, which also verified the Re_R2007 scheme's advantages. Finally, in comparison with other global advection schemes, the Re_R2007 scheme was competitive in terms of accuracy and flow independence. An encouraging possibility for the application of the Re_R2007 scheme to the GRAPES model is provided.

摘要: The Global/Regional Assimilation and PrEdiction System (GRAPES) is the new-generation numerical weather prediction (NWP) system developed by the China Meteorological Administration. It is a fully compressible non-hydrostatical global/regional unified model that uses a traditional semi-Lagrangian advection scheme with cubic Lagrangian interpolation (referred to as the SL_CL scheme). The SL_CL scheme has been used in many operational NWP models, but there are still some deficiencies, such as the damping effects due to the interpolation and the relatively low accuracy. Based on \hboxReich's semi-Lagrangian advection scheme (referred to as the R2007 scheme), the Re_R2007 scheme that uses the low- and high-order B-spline function for interpolation at the departure point, is developed in this paper. One- and two-dimensional idealized tests in the rectangular coordinate system with uniform grid cells were conducted to compare the Re_R2007 scheme and the SL_CL scheme. The numerical results showed that: (2) the damping effects were remarkably reduced with the Re_R2007 scheme; and (3) the normalized errors of the Re_R2007 scheme were about 7.5 and 3 times smaller than those of the SL_CL scheme in one- and two-dimensional tests, respectively, indicating the higher accuracy of the Re_R2007 scheme. Furthermore, two solid-body rotation tests were conducted in the latitude-longitude spherical coordinate system with non-uniform grid cells, which also verified the Re_R2007 scheme's advantages. Finally, in comparison with other global advection schemes, the Re_R2007 scheme was competitive in terms of accuracy and flow independence. An encouraging possibility for the application of the Re_R2007 scheme to the GRAPES model is provided.

1. Introduction
  • The transport process is of crucial importance in atmospheric movements. In numerical weather prediction (NWP) models, some discretized advection schemes need to be used to approximate the solution of the advection equations. The accuracy of the advection scheme plays a key role in the performance of NWP models. Therefore, lots of attention has been paid to the study of advection schemes in the past decades. In the early NWP models, the Eulerian advection scheme was widely used, but its time step length is strictly limited by the Courant-Friedrichs-Lewy condition for ensuring the numerical stability of integration in NWP models. In recent years, semi-Lagrangian advection schemes incorporating the semi-implicit scheme have been widely used in operational weather and climate models because they use a larger time step length than Eulerian advection schemes, which reduces the time taken to run the model and lowers the computational cost. Robert (1981, 1982) proved that the combination of the semi-Lagrangian treatment of advection and the semi-implicit treatment of gravitational oscillations could increase the maximum stable time step length by a factor of six at some extra cost but without lowering the numerical accuracy. In addition,(Staniforth and C#cod#x0020f;t#cod#x000e9;, 1991) demonstrated that the semi-Lagrangian advection scheme could achieve a computational accuracy comparable to the Eulerian scheme and that it also exhibited better efficiency and less dissipation of quantity. As expected, these advantageous semi-Lagrangian advection schemes have been introduced into many NWP models at operational meteorological centers, such as the European Center for Medium-range Weather Forecasts (ECMWF), the UK Met Office, the METEO FRANCE, the National Centers for Environmental Prediction (NCEP) in USA, Canadian Meteorological Center (CMC), and Japan Meteorological Agency (JMA), etc.

    GRAPES (Global/Regional Assimilation and PrEdiction System), the new-generation NWP model developed by the China Meteorological Administration (CMA), adopted the traditional semi-Lagrangian advection scheme with cubic Lagrangian interpolation (SL_CL scheme), since the cubic Lagrangian interpolation represents a good compromise between accuracy and computational cost (Staniforth and C#cod#x0020f;t#cod#x000e9;, 1991). However, the SL_CL scheme damps the extrema of quantities and incurs some spurious oscillations in the discontinuous areas due to the interpolation (McDonald, 1984, 1987), so it cannot preserve the shape of the initial field. This behavior also indicates a potential problem in the traditional semi-Lagrangian advection scheme. The simple interpolation (e.g., linear interpolation) can reduce the computational cost and smooth the spurious oscillations of quantities in the discontinuous areas but sacrifices the numerical accuracy; while the complex interpolation (e.g., the high-order interpolation method) has higher accuracy but with expensive computational costs and severe spurious oscillations of quantities in the discontinuous areas. To suppress the spurious oscillations in the discontinuous areas, (Bermejo and Staniforth, 1992) converted the traditional semi-Lagrangian advection scheme into the quasi-monotone semi-Lagrangian (QMSL) scheme by combining the low- and high-order approximate solutions. This kind of treatment in the QMSL scheme can effectively control the oscillations and be easily incorporated into any traditional semi-Lagrangian advection scheme. However, it should be noted that it reduces the numerical accuracy in the QMSL scheme. Therefore, it is expected to improve the numerical accuracy and simultaneously reduce the spurious oscillations in the discontinuous areas in the semi-Lagrangian advection scheme.

    Based on the remapped particle-mesh (RPM) method and the forward trajectory algorithm, (Cotter et al., 2007) developed the remapped particle-mesh semi-Lagrangian (RPM SL) advection scheme from the continuity equation, in which the factor of grid cell area was introduced to the interpolation function to ensure the conservativeness property of this scheme. It was found that the RPM SL scheme could theoretically keep the conservativeness property of quantities and simultaneously achieve the same accuracy as the SL_CL scheme. However, for the RPM SL scheme, a tri-diagonal linear system of equations needs to be solved at each time step, and this increases the computational cost, especially when the resolution is increased. To overcome this deficiency, (Reich, 2007) employed a quasi-interpolation (Powell, 1981) to approximate the solution of the tri-diagonal linear system of equations in the RPM SL scheme and advanced the explicit RPM SL scheme based on the forward trajectory algorithm, which, as has been proved, reduces the computational cost without sacrificing the conservativeness property and numerical accuracy of the RPM SL scheme. At the end of Reich's paper (see section 6 of his paper), he appropriately adjusted the interpolation algorithm in the explicit RPM SL scheme, and applied it to the backward trajectory algorithm (hereafter, referred to as the R2007 scheme), but he did not verify the feasibility of the R2007 scheme in practical tests.

    During our tests of the R2007 scheme, it was found that on the grid mesh with uniform grid cells, the factor of grid cell area had nothing to do with the whole computation since it was offset in the interpolation process. Also, when the R2007 scheme was applied to the grid mesh with non-uniform grid cells (e.g., the latitude-longitude spherical grid mesh), this factor of grid cell area could not guarantee the conservativeness property in theory but also possibly led to numerical instability in the interpolation (e.g., when the fluid flows across the poles). In this paper, we attempt to remove the factor of grid cell area in the R2007 scheme and only consider the effect of the grid spacing in the interpolation functions to make the R2007 scheme valid in a latitude-longitude grid coordinate system. For convenience, the R2007 scheme after this revision is referred to as the Re_R2007 scheme. The feasibility of the Re_R2007 scheme was justified and compared with the SL_CL scheme on grid meshes with uniform and non-uniform grid cells and the results were promising and encouraging.

    The paper is organized as follows: Section 2 briefly discusses the semi-Lagrangian method and cubic Lagrangian interpolation. Section 3 is devoted to the description and analysis of the specific procedures of the Re_R2007 scheme. In section 4, the results from one-dimensional (1D) and two-dimensional (2D) numerical tests in the rectangular coordinate system with uniform grid cells are presented and analyzed to compare the Re_R2007 and SL_CL schemes. Section 5 gives the results of solid-body rotation tests in the latitude-longitude spherical coordinate system with non-uniform grid cells in both schemes. Finally, some discussions and conclusions about the Re_R2007 scheme are drawn in section 6.

2. Semi-Lagrangian method and cubic Lagrangian interpolation
  • Detailed descriptions of GRAPES can be seen in the work of (Xue, 2004), (Xue and Liu, 2007), (Chen et al., 2008), (Yang and Shen, 2011) and (Liu et al., 2012). The key project for the development of GRAPES was launched in 2001, and it has been operationally or quasi-operationally implemented since 2004 at the national and regional meteorological centers of the CMA. GRAPES uses the SL_CL scheme to deal with the solution of the advection equations. For convenience, the discussion in this section focuses on the solution of the pure advection problem (i.e., the forcing term in the advection equation was ignored) in the 2D case. Then, the Lagrangian advection equation for the scalar f(x,t) can be written as

    where denotes the position vector; t the time; V the wind field; and #cod#8711; the gradient operator; and d/dt and #cod#8706;/#cod#8706;t are the Lagrangian and Eulerian time derivatives, respectively.

    The semi-Lagrangian approximation to Eq. (2) is generally expressed as

    where t(n)=n#cod#916; t and t(n+1)=(n+1)#cod#916; t; n represents the number of the time step; and #cod#916; t is the time-step length. The subscripts "D" and "A" denote the departure point at time t(n) and the arrival point at time t(n+1), respectively.

    Equation (4) reveals that for the pure advection problem, the value at the arrival point at time t(n+1) is equal to that at the departure point at time t(n). When it is applied on a regular grid mesh, the grid point is assumed to be the arrival point at time t(n+1). Before the value at the grid point at time t(n+1) can be computed, the position of its corresponding departure point at time t(n) must be located and its value estimated. This is the essence of the semi-Lagrangian method. In general, the departure point does not match a regular grid point; its value needs to be computed by some interpolation method. In GRAPES, the cubic Lagrangian interpolation method is used and its interpolation procedures are described below.

    The 1D cubic Lagrangian interpolation polynomial is expressed as

    f(x)= l0(x)f0+l1(x)f1+l2(x)f2+l3(x)f3, (4)

    where x denotes the interpolated point; xi(i=0,1,2,3) are the four distinct known points surrounding the interpolated point; and fi=f(xi) and li(x) represent the value and cubic Lagrangian basis function at the corresponding point. Accordingly, li(x) can be written as

    In the 2D case, 16 points around the interpolated point are required to compute the value at the interpolated point, according to (Li et al., 2006). Figure 1 shows a regular grid mesh with uniform grid spacing that was introduced to illustrate the specific procedures of the 2D cubic Lagrangian interpolation method, where the Eulerian grid point and its corresponding value are denoted by xi,j=(xi,yj)=(i#cod#183;#cod#916; x,j#cod#183;#cod#916; y) and fi,j=f(xi,yj),i=1,#cod#8230; M,j=1,#cod#8230;,N,M(N) and #cod#916; x(#cod#916; y) are the total grid number and the uniform grid spacing in x(y) direction. Point A represents any Eulerian grid point and signifies the arrival point in the semi-Lagrangian method, and point D is thought to be its corresponding departure point denoted by x D=(x D,y D). The value of f at point D can be interpolated from the values of f at the surrounding 16 points enclosed by the thick solid lines in Fig. 1. First, the values of f at points a, b, c and d are computed using the 1D cubic Lagrangian interpolation along the x direction, then the value of f at point D is determined using the 1D cubic Lagrangian integration with the values of f at points a, b, c and d. The detailed interpolation procedures are presented below:

    f(xa,ya) = la0fi−1,j−1+la1fi,j−1+la2fi+1,j−1+la3fi+2,j−1, (9)

    f(xb,yb) = lb0fi−1,j+lb1fi,j+lb2fi+1,j+lb3fi+2,j, (10)

    f(xc,yc) = lc0fi−1,j+1+lc1fi,j+1+lc2fi+1,j+1+lc3fi+2,j+1, (11)

    f(xd,yd) = ld0fi−1,j+2+ld1fi,j+2+ld2fi+1,j+2+ld3fi+2,j+2, (12)

    f(xD,yD) = lD0f(xa,ya)+lD1f(xb,yb)+lD2f(xc,yc)+lD3f(xd,yd), (13)

    where lai,lbi,lci,l di and l Di(i=0,1,2,3) are the cubic Lagrangian basis functions corresponding to the points a, b, c, d and D.

    Figure 1.  Interpolation schematic for the SL_CL and Re_R2007 schemes.

3. Interpolation algorithm in the Re_R2007 scheme
  • The Re_R2007 scheme was developed during tests of the R2007 scheme. It removes the factor of the grid cell area in the whole computation in the R2007 scheme and only allows for the effect of the grid spacing in the basis functions. So its main approach and general procedure are consistent with those of the R2007 scheme and can be seen in the work of (Reich, 2007). According to the work of (Cotter et al., 2007) and (Reich, 2007), in this section we briefly describe the interpolation algorithm of the Re_R2007 scheme in the 2D case in the context of Fig. 1. It can be easily modified to the 1D form by adjusting the dimensions in the computation. In fact, the Re_R2007 scheme associates the cubic B-spline function with the linear B-spline function and takes advantage of the low-order value to correct the high-order value, contributing to the higher accuracy of the Re_R2007 scheme. In Fig. 1, each Eulerian grid point xi,j, carries a bi-cubic and bi-linear B-spline function, which are defined as, respectively,

    where x=(x,y) denotes the interpolated point and #cod#968; cs(r) and #cod#968; ls(r) are the cubic and linear B-spline basis function, respectively:

    In Fig. 1, M#cod#215; N Lagrangian particles are introduced in the Re_R2007 scheme. Based on the remapped particle-mesh (RPM) method (Cotter et al., 2007), the M#cod#215; N Lagrangian particles are uniformly assigned onto the M#cod#215; N Eulerian grid points at the beginning of each time step. In the meantime, the value of the Lagrangian particles is defined as Fi,j=F(xi, yj), which is unknown.

    To estimate the value of F of the Lagrangian particles, it is assumed that the value of f at the Eulerian grid points and the value of F of the Lagrangian particles are related in the following way:

    Equation (19) indicates that the value of f at any Eulerian grid point can be interpolated from the values of F at the surrounding 9 Lagrangian particles. In this case, considering the values of f at all the Eulerian grid points (the exact value at the Eulerian grid points, referred to as f ex), Eq. (19) can constitute a linear system of equations about the value of F. By solving this linear system of equations, we can obtain the value of F of the Lagrangian particles (the exact value of the Lagrangian particles, referred to as F ex). Unfortunately, this treatment causes expensive computational cost. To avoid this, a quasi-interpolation (Powell, 1981) is employed to approximate the value of F of the Lagrangian particles (the approximate value of the Lagrangian particles, referred to as F ap), which is expressed as

    However, it should be noted that there exists an error between Fex from Eq. (19) and Fap from Eq. (20), so when F ap is put into Eq. (19) to compute the value of f at the Eulerian grid points (the approximate value at the Eulerian grid points, referred to as f ap), there is an error between f ex and fap at the Eulerian grid points, which is defined as #cod#916; f at any Eulerian grid point:

    To eliminate the effects of the error between F ex and F ap, the linear combination of Eqs. (20) and (21) are used to compute the value of f at point D:

    Eqs. (20)-(22) describe the interpolation procedures in the Re_R2007 scheme, the final value of f at point D can be divided into two parts: one part is the first term on the right-hand side of Eq. (22), which is computed using the bi-cubic B-spline function with the values of F at the surrounding 16 Lagrangian particles (i.e., the 16 points enclosed by the thick solid lines in Fig. 1); the other part is the second term on the right-hand side of Eq. (22), which is interpolated using the bi-linear B-spline function with the values of #cod#916; f at the surrounding four Eulerian grid points (i.e., the four points enclosed by the thin solid lines in Fig. 1). The numerical tests showed that, in practice, the second part is mainly used to correct the first part, thus contributing to the improvement in the accuracy of the Re_R2007 scheme. When the field of quantities is sufficiently smooth, the correction of the second part to the first part is very small, since the value of Fap from Eq. (20) closely approximates the value of Fex from Eq. (19), then the value of #cod#916;f can be almost negligible so that the high-order interpolated value is sufficient to approximate the value of f at point D. However, if the field of quantities is not sufficiently smooth, then the second part can influence the first part positively and significantly, since the value of #cod#916; f can be very important due to the relatively big error between the value of F ex from Eq. (19) and the value of F ap from Eq. (20), then the linear interpolation with the value of #cod#916;f can make a big difference in correcting the high-order interpolated results.

    By comparing the interpolation procedures between the SL_CL scheme and the Re_R2007 scheme, it was found that the Re_R2007 scheme increased the computational cost relatively for the extra efforts to compute F and #cod#916;f before the interpolation, while in the SL_CL scheme the value at the departure point could be directly interpolated from the values of f at the surrounding grid points. However, the numerical tests showed that the extra computational cost in the Re_R2007 scheme only increased a little more and could be acceptable under the existing computational capability. On the other hand, its remarkable accuracy would overcome the deficiency in this aspect. Therefore, the Re_R2007 scheme should be a good scheme to be incorporated into the numerical model.

4. Numerical tests in the rectangular coordinate system with uniform grid cells
  • To evaluate the behavior of the Re_R2007 scheme, we first compared it with the SL_CL scheme using two groups of numerical tests: the 1D and 2D idealized tests in the rectangular coordinate system with uniform grid cells.

    For the sake of contrasting the difference between the numerical solution and the exact solution, we defined the numerical error as #cod#963;:

    #cod#963;i,j=(fi,j)num−(fi,j)exa, (22)

    where the subscripts (i,j) represents the ith and jth Eulerian grid point, and the subscripts "num" and "exa" denote the corresponding numerical solution and exact solution, respectively.

    Meanwhile, to better illustrate the numerical accuracy, the normalized errors were defined according to (Williamson et al., 1992):

    where the relative parameters are the same as Eq. (23) and I(f) represents the domain integral defined by

    where M and N are the total number of grid points in the x and y direction, respectively.

  • Considering that the interpolation can incur some damping to the extrema of quantities, a sine wave was chosen in the 1D test to examine the damping effects in both schemes by inspecting the effects on the peak and valley of a sine wave.

    The initial distribution was:

    f(x,0)=sin(#cod#x003c0;x), −1.0#cod#x02264;x#cod#x02264;1.0 . (27)

    In the test, the uniform grid spacing was #cod#916; x=0.02, the time step was #cod#916; t=0.01, and the constant velocity was u=1.0. The number of time steps was 2000 with the completion of 10 revolutions.

    Figure 2 illustrates the numerical results (top figures) and errors (bottom figures) of the SL_CL and Re_R2007 schemes. It was found that both advection schemes could favorably maintain the initial distribution of the sine wave with negligible deformation after a 2000-step integration. However, the peak and valley of the sine wave in the SL_CL scheme were more seriously damped than in the Re_R2007 scheme. Figs. 2c and d show that the high-value errors mainly rest near the peak and valley of the sine wave, and the maximum of the absolute errors in the SL_CL scheme is roughly 0.0007, while it is not more than 0.0001 in the Re_R2007 scheme.

    Figure 2.  Numerical results and errors of the 1D sine wave experiment. Numerical results of the (a) SL_CL scheme and (b) Re_R2007 scheme; and numerical errors of the (c) SL_CL scheme and (d) Re_R2007 scheme.

    Figure 3 illustrates the normalized errors of both schemes, showing the distinctly higher accuracy of the Re_R2007 scheme. In both schemes, the three normalized errors increased at nearly the same rate, but the increase rate in the Re_R2007 scheme was much smaller than that in the SL_CL scheme. At the end of the integration, the ratio of the normalized errors in the SL_CL scheme to those in the Re_R2007 scheme was approximately 0.00075/0.0001=7.5. Therefore, the Re_R2007 scheme showed a 7.5-time higher accuracy than the SL_CL scheme in this respect.

    Figure 3.  Normalized errors of the 1D sine wave experiment: (a) SL_CL scheme; (b) Re_R2007 scheme (solid line: L1; dotted line: L2; dashed line: L#cod#8734;).

  • The interpolation also easily incurs spurious oscillations in the discontinuous areas, affecting the numerical accuracy. For instance, in the initial field with all positive values, negative values may occur after the interpolation. In the 2D test, we chose a cosine bell with an unsmooth boundary and all non-negative values, to examine the spurious oscillations as well as the damping effects of both schemes in the 2D case.

    The initial distribution of the cosine bell is given by

    where .

    The test settings were as follows: the uniform grid spacing #cod#916; x=#cod#916; y=0.02; the integration domain [-1.0,1.0]#cod#215;[-1.0,1.0], the time step #cod#916; t=0.01; and the constant 2D velocity u=v=1.0. This test was run with 2000 time steps to complete 10 revolutions.

    The numerical results (see Figs. 4a and b) indicate that both schemes developed obvious spurious oscillations near the unsmooth boundary. However, the size of the oscillation in the Re_R2007 scheme was smaller than that in the SL_CL scheme. This was also seen in the corresponding numerical error figures (see Figs. 4c and d). In the SL_CL scheme, there were several high-value areas of positive/negative errors ranging from -0.016 to 0.016 near the unsmooth boundary and the errors inside the cosine bell (namely, the smooth areas) were also very serious; while the Re_R2007 scheme only had small errors ranging from -0.004 to 0.004 near the unsmooth boundary and the errors inside the cosine bell were negligible. Hence, the Re_R2007 scheme showed marked advantages over the SL_CL scheme in reducing the spurious oscillations in the unsmooth boundary.

    Figure 4.  Numerical results and errors of the 2D cosine bell experiment. Numerical results of the (a) SL_CL scheme and (b) Re_R2007 scheme; and numerical errors of the (c) SL_CL scheme and (d) Re_R2007 scheme.

    To further inspect the damping effects in both schemes, we compared the cross sections at y=0 of the numerical solution and the exact solution. Here, we defined the cross section errors as the difference of the cross sections between the numerical solution and the exact solution. As is clearly shown in Fig. 5, negative values occur near the unsmooth boundary of the cosine bell in both schemes. Meanwhile, it should be noted that the maximum of the cosine bell in the SL_CL scheme was not more than 0.99 after a 2000-step integration, while the effect on this maximum in the Re_R2007 scheme was negligible. Hence, it is easily concluded that the damping effects in the Re_R2007 scheme are dramatically weakened in the 2D case compared with the SL_CL scheme.

    Figure 5.  Cross section and errors at Y = 0 of the 2D cosine bell experiment. Cross section of the (a) SL_CL scheme and (b) Re_R2007 scheme (solid line: cross section of the numerical results; dotted line: cross section of the initial field). Cross section errors of the (c) SL_CL scheme and (d) Re_R2007 scheme.

    Figure 6 shows the normalized errors of both schemes. Unlike the 1D test, L1,L2 and L#cod#8734; in the 2D test increased at different rates along the integration time, and they increased quickly at first but gently afterwards. The increase rate of the normalized errors in the Re_R2007 scheme was consistently much smaller than those in the SL_CL scheme. For example, at the end of the integration L1,L2 and L#cod#8734; in the SL_CL scheme were approximately 0.038, 0.023 and 0.016, respectively, while the corresponding values in the Re_R2007 scheme were roughly 0.008, 0.006 and 0.0045. The normalized errors in the SL_CL scheme were more than three times larger than those in the Re_R2007 scheme. Therefore, the Re_R2007 scheme achieved a distinctly higher accuracy.

    Figure 6.  Normalized errors of the 2D cosine bell experiment: (a) SL_CL scheme and (b) Re_R2007 scheme (solid line: L1; dotted line: L2; dashed line: L#cod#8734;).

    From the 1D and 2D idealized tests in the rectangular system with uniform grid cells, it can be concluded that the Re_R2007 scheme advantageously reduces the spurious oscillations of quantities and the damping effects due to the interpolation as well as significantly improves the numerical accuracy.

5. Solid-body rotation tests in the spherical coordinate system with nonuniform grid cells
  • A 2D advection test of a solid body with a divergence-free current was suggested by (Williamson et al., 1992), and it is now commonly used to evaluate global advection schemes. In this section, we report the results from performing solid-body rotation tests in the latitude-longitude spherical coordinate system with non-uniform grid cells, to evaluate the behavior of the Re_R2007 scheme. For comparison, we also make use of such terms as the numerical errors and normalized errors defined by Eqs. (23)-(26), in which the domain integral I(f) is replaced by the global integral

    where #cod#923; and #cod#966; denote the longitude and latitude, respectively.

    The initial distribution was a cosine bell like the 2D idealized test, and in the latitude-longitude coordinate system it can be written as (see Fig. 7):

    where R=1/3, and r is the distance between (#cod#923;,#cod#966;) on a great circle, and the center of the cosine bell is (#cod#923;0,#cod#966;0)=(#cod#960;/2,0).

    The velocity components of the advecting wind field are given by

    u=u0(cos#cod#966;cos#cod#945;+sin#cod#966;cos#cod#923;sin#cod#945;), (32)

    v=−u0sin#cod#923;sin#cod#945;, (33)

    where u0=2#cod#960;/(12 days), and #cod#945; is the angle between the axis of the solid-body rotation and the polar axis of the spherical coordinate system. The tests were run at #cod#945;=0 and #cod#945;=#cod#960;/2, i.e., the cosine bell traveled along the Equator and across the poles, respectively.

    The spherical (#cod#923;,#cod#966;) domain consisted of a 128#cod#215; 64 uniform-resolution (2.8125#cod#x000b0;) mesh. The time step (4050 s) was chosen such that the solid body took 256 time steps to complete one revolution around the globe. To accurately locate the position of the departure point, the procedures of the Ritchie and Beaudoin algorithm (Ritchie and Beaudoin, 1994) were used in the low-latitude regions (between 80#cod#x000b0;S and 80#cod#x000b0;N), while the method of McDonald and Bates (1989) was employed in the high regions (over 80#cod#x000b0;S and 80#cod#x000b0;N) where the Ritchie and Beaudoin algorithm breaks down due to the serious curvature effects near the poles.

    Figure 7.  Initial distribution of the solid-body rotation tests.

    Figures 8 and 9 show the numerical results and errors of both schemes at #cod#945;=0 and #cod#945;=#cod#960;/2. It is noted that the cosine bell has undergone a distinct stretching in the flow direction in the SL_CL scheme at #cod#945;=0 and #cod#945;=#cod#960;/2, while the stretching effect in the Re_R2007 scheme is not so serious. Meanwhile, the maximum value of the cosine bell in the SL_CL scheme was damped to 0.8 after one revolution, but this value in the Re_R2007 scheme remained at 0.9 like the initial field. This distinction was also seen in the values of the numerical errors (see Fig. 9). The errors mainly occurred along the flow direction with several high-value areas of positive/negative errors, so we can observe the spurious oscillations near the unsmooth boundary along the flow direction. For the same scheme (e.g., SL_CL or Re_R2007 scheme), the numerical errors were generally at the same level whenever the cosine bell traveled along the equator or across the poles. For example, the numerical errors range from -0.1 to 0.1 in the SL_CL scheme, while they range from -0.02 to 0.02 in the Re_R2007 scheme. In addition, for the flow in the same direction (e.g., #cod#945;=0 or #cod#945;=#cod#960;/2), the intensity and size of the numerical errors in the Re_R2007 scheme are rather smaller than those in the SL_CL scheme. For instance, the maximum of the absolute errors of the Re_R2007 scheme (roughly 0.02) is only about one-fifth of that of the SL_CL scheme (roughly 0.1). Therefore, on the spherical grid mesh with non-uniform grid cells, the Re_R2007 scheme can also effectively reduce the spurious oscillations in the discontinuous areas and the effect of damping due to the interpolation.

    Figure 8.  Numerical results of the solid-body rotation tests (#cod#945;=0 and #cod#945;=#cod#960;/2): #cod#945;=0 for the (a) SL_CL scheme and (b) Re_R2007 scheme; #cod#945;=#cod#960;/2 for the (c) SL_CL scheme and (d) Re_R2007 scheme.

    Figure 9.  Numerical errors of the solid-body body rotation tests (#cod#945;=0 and #cod#945;=#cod#960;/2): #cod#945;=0 for the (a) SL_CL scheme and (b) Re_R2007 scheme; #cod#945;=#cod#960;/2 for the (c) SL_CL scheme and (d) Re_R2007 scheme.

    Figure 10.  Normalized errors of the solid-body rotation tests: #cod#945;=0 for the (a) SL_CL scheme and (b) Re_R2007 scheme; #cod#945;=#cod#960;/2 for the (c) SL_CL scheme and (d) Re_R2007 scheme (solid line: L1; dotted line: L2; dashed line: L#cod#8734;).

    Figure 10 shows the normalized errors of both schemes at #cod#945;=0 and #cod#945;=#cod#960;/2. These results are similar to those in the 2D numerical test, and the normalized errors of the SL_CL scheme are approximately three times larger than those of the Re_R2007 scheme at #cod#945;=0 and #cod#945;=#cod#960;/2, so the Re_R2007 scheme has a distinctly higher accuracy. In addition, to further verify the Re_R2007 scheme, we compared our numerical results from the solid-body rotation tests with those of other existing global advection schemes, such as the semi-Lagrangian transport on a reduced grid (RG, Rasch, 1994), the flux-form semi-Lagrangian scheme (FFSL, Lin and Rood, 1996), the cell-integrated semi-Lagrangian scheme (CISL, Nair and Machenhauer, 2002) and the semi-Lagrangian inherently conserving and efficient scheme (SLICE, Zerroukat et al., 2004). For convenience, we directly referred to the numerical results of those schemes [see Table 1 and 2 of (Nair and Machenhauer, 2002), Table 1 of (Zerroukat et al., 2004), and Table 4 of (Li et al., 2008)], in which the MAX and MIN errors are defined as (Rasch, 1994). Table 1 shows that the Re_R2007 scheme presents significant advantages over the SL_CL and RG schemes in all the error terms except the MIN error of RG2.8-m at #cod#945;=#cod#960;/2. Compared with the FFSL-3, FFSL-5 and CISL-M schemes, the Re_R2007 scheme also shows definite superiorities with respect to the L2, L#cod#8734; and MAX terms (expect for the L2 term of FFSL-5 at #cod#945;=#cod#960;/2). On the contrary, the Re_R2007 scheme cannot effectively compete with the CISL-P and SLICE-S schemes in most of the error terms. However, the accuracy of the Re_R2007 scheme is hardly affected by the flow directions, which is indicated by the fact that the Re_R2007 scheme can achieve nearly the same level of the normalized errors at #cod#945;=0 and #cod#945;=#cod#960;/2, while this is not the case for the CISL-P and SLICE-S schemes.

6. Discussion and conclusions
  • In this paper, the Re_R2007 scheme with the low- and high-order B-spline interpolation function, in which a correction term is added by the low-order interpolation and has a significant effect on the improvement of its numerical accuracy, has been presented. To evaluate the Re_R2007 scheme effectively, comparisons with the SL_CL scheme in GRAPES were made using 1D and 2D idealized tests in the rectangular coordinate system with uniform grid cells, as well as two solid-body rotation tests in the latitude-longitude spherical coordinate system with non-uniform grid cells. It was demonstrated that the Re_R2007 scheme shows significant advantages over the SL_CL scheme and the conclusions are as follows:

    (2) The Re_R2007 scheme can effectively reduce the spurious oscillations near the discontinuous areas as well as the damping effects due to the interpolation so that it can better maintain the initial distribution of quantities. In all the tests, the intensity and size of the spurious oscillations near the discontinuous areas in the Re_R2007 scheme were dramatically smaller than the SL_CL scheme. Furthermore, the extrema were damped a little or even negligibly in the Re_R2007 scheme while the damping effects in the SL_CL scheme were relatively serious.

    (3) The Re_R2007 scheme can achieve an obviously higher accuracy in comparison with the SL_CL scheme. At the end of the integration, the ratios of the normalized errors in the SL_CL scheme to those in the Re_R2007 scheme were roughly 7.5 and 3 in the 1D and 2D idealized tests, respectively, indicating that the Re_R2007 scheme has a higher accuracy.

    (4) In the solid-body rotation tests in the latitude-longitude spherical coordinate system with non-uniform grid cells, the Re_R2007 scheme exhibited similar advantages over the SL_CL scheme as in the 2D idealized test, achieving a nearly three-times higher accuracy. In addition, by comparing the Re_R2007 scheme with other existing global advection schemes, it can be seen that the Re_R2007 scheme can be competitive with some advection schemes in terms of accuracy. It should be noted, at the same time, that the accuracy of the Re_R2007 scheme was hardly affected by flow directions in the tests, but this is not the case for the CISL-P and SLICE-S schemes.

    In the near future, the Re_R2007 scheme will be incorporated into the operational NWP system of GRAPES at the China Meteorological Administration. It is expected that some 3D idealized tests and a series of simulations of real cases will be performed to further investigate the feasibility of the Re_R2007 scheme in GRAPES.

Reference

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint