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.