Quarterly Journal of the Royal Meteorological Society

Q. J. R. Meteorol. Soc. 137: 1913–1926, October 2011 A

The Canadian Global Environmental Multiscale model on the Yin-Yang grid system Abdessamad Qaddouri* and Vivian Lee Meteorological Research Division, Atmospheric Science and Technology Directorate, Environment Canada, Dorval, Qu´ebec, Canada *Correspondence to: A. Qaddouri, Environment Canada, 2121 Trans-Canada Hwy, Dorval QC H9P 1J3, Canada. E-mail: [email protected]

At the Canadian Meteorological Center (CMC), we are currently developing the future global forecasting Yin-Yang model. In the horizontal we use spherical coordinates on the overset Yin-Yang grid, while in the vertical we use a loghydrostatic-pressure coordinate on the Charney–Phillips grid. The parametrization of physical processes is kept the same as in the current Global Environmental Multiscale (GEM) operational model. The Yin-Yang global forecast is performed by considering a domain decomposition (a two-way coupling method) between two limited-area models (LAMs) discretized on the two panels of the Yin-Yang grid and using the same time step. Each panel of the Yin-Yang grid system is extended by a static halo region and uses the same fully implicit semi-Lagrangian method as in the GEM operational model to solve its own dynamic core. The spatial and time discretizations are implemented independently on each quasi-uniform latitude–longitude subgrid. The static halo region plays the same role as the piloting region in limited-area modelling. Since the two subgrids of the Yin-Yang grid do not match, the update of the variables in the pilot region is done by cubic Lagrange interpolation. For our model validation, we ran 42 winter and 42 summer cases using analysis from 2008–2009 and we compared five-day forecast results against observations. No noise is seen in the overlap regions during the simulations. Preliminary results presented in this article are encouraging and demonstrate that in comparison with observations the new Yin-Yang system performs as well as the GEM global c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd. model. Key Words: method

overset grid on the sphere; semi-Lagrangian semi-implicit scheme; domain decomposition

Received 22 December 2010; Revised 3 May 2011; Accepted 25 May 2011; Published online in Wiley Online Library 1 August 2011 Citation: Qaddouri A, Lee V. 2011. The Canadian Global Environmental Multiscale model on the Yin-Yang grid system. Q. J. R. Meteorol. Soc. 137: 1913–1926. DOI:10.1002/qj.873

1. Introduction The deterministic Global Environmental Multiscale (GEM) operational model (Yeh et al., 2002) constitutes a very important tool for medium-range numerical weather prediction at the Canadian Meteorological Center (CMC). Over the last decade, significant improvements have been made in the GEM physical process representation (B´elair

et al., 2009) and the model lid has been raised to 0.1 mb with no major changes in its dynamical core part. In the current year, however, the non-staggered vertical grid with a hydrostatic pressure coordinate in the GEM model (Yeh et al., 2002) was replaced by the Charney–Phillips grid with a log-hydrostatic-pressure coordinate (Girard et al., 2010). This new spatial vertical discretization has a positive impact on model performance, since it improves the accuracy of

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

1914

A. Qaddouri and V. Lee

area modelling (LAM). Each subdomain is extended by a static halo region that plays the same role as the piloting region in the LAM. The Schwarz method consists of solving the subdomain problems iteratively by using the subdomain solutions to update the boundary conditions in the pilot region of the neighbouring subdomain. Since the two subgrids of the Yin-Yang grid do not match, a full cubic Lagrange interpolation is used to bring the boundary conditions to each respective target grid. This technique has been used to solve the nonlinear shallow-water equations on the Yin-Yang grid in Qaddouri (2011). The three-dimensional (3D) dynamical core has already been developed and tested on both the polar-cap grid (Dudhia and Bresch, 2002) and the cubed sphere grid (Adcroft et al., 2004). The Yin-Yang grid has just been considered for the first time in the development of a 3D dynamical core for simulating analytical cases of atmospheric flow ( Baba et al., 2010). In this article, we will focus on the feasibility of the Yin-Yang grid for operational medium-range global forecasts while considering both dynamical and physical process representations. The model formulation and numerical solution are described in section 2. For simplicity, only hydrostatic model formulations will be shown, as the non-hydrostatic equations involve a few extra terms but do not change the technique used in obtaining the numerical solution. The schemes used for the parametrization of physical processes are presented in section 3. The use of a digital filter and horizontal diffusion is described in section 4. The numerical results are shown and discussed in section 5. Finally, our conclusions are presented in section 6.

(a)

(b)

Figure 1. (a) Latitude–longitude grid and (b) Yin-Yang grid. This figure is available in colour online at wileyonlinelibrary.com/journal/qj

2. 2.1.

the hydrostatic and non-hydrostatic equations by defining the vertical finite differences at logarithmic midpoints. This staggered vertical grid GEM version is expected to become operational in 2011, and we will refer to it as the GEM global model. A problem associated with the construction of any global atmospheric model is the numerical representation of the spherical geometry. Currently, the GEM global model uses the latitude–longitude (denoted Lat–Lon) grid system, which leads to singularities and convergence of meridians in the polar regions. Using today’s computer architectures with distributed memory, this may even result in an unbalanced computational load in the context of domain decomposition with message-passing interfaces. In this paper, we developed a Yin-Yang model by replacing a Lat–Lon mesh grid as in Figure 1(a) with a Yin-Yang grid system (Figure 1(b)) that imposes a more uniform grid spacing. The Yin-Yang model uses the Schwarz domain decomposition method (DDM) in order to solve a system of forced primitive equations (PEs) on a quasi-uniform overset grid free from polar singularity. This grid system is constructed by overlapping two perpendicularly oriented identical parts of a global Lat–Lon grid (Figure 1(b)) as suggested by Kageyama and Sato (2004). In our implementation of the DDM, the governing equations on each (Yin or Yang) subgrid are the PEs and these are discretized using the implicit semi-Lagrangian method as in the GEM operational model (Cˆot´e et al., 1998; Yeh et al., 2002). The solution for each subdomain is based on the same technique used in limited c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Model formulation Governing equations

In each of the two subdomains of the Yin-Yang grid system, the governing equations are the forced hydrostatic primitive equations described in Girard et al. (2010) as dVh + f k × Vh + RT∇ζ Bs + ∇ζ φ = FH , dt T d d ˙ ln −κ (Bs) + ζ = F T , dt T∗ dt ∂B d ∂ ζ˙ Bs + ln 1 + s + ∇ζ · Vh + + ζ˙ = 0, dt ∂ζ ∂ζ ∗ ∂ ζ − φ /RT T = 0, − T∗ ∂ (ζ + Bs)

(1) (2) (3)

(4)

where ζ = ln p − Bs is a log-hydrostatic-pressure type coordinate with p the pressure, B = [(ζ − ζtop )/(ζsurf − ζtop )] the metric term, s = ln psurf − ζsurf the mass variable that depends only on the horizontal, T the temperature, φ = φ ∗ + φ the geopotential, φ ∗ (ζ ) = −RT ∗ (ζ − ζsurf ), ∗ T = const and Vh the horizontal wind. R is the gas constant and κ = R/cp , where cp is the specific heat at constant pressure. The coordinate system permits a switch-controlled choice between the hydrostatic primitive equations and non-hydrostatic Euler equations. See Girard et al. (2010) for additional details. The prognostic equations (1)–(3) are respectively the horizontal momentum, thermodynamic and continuity Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

equations, and (4) is the diagnostic hydrostatic equation. The terms FH and F T are parametrized physical forcings. It is also convenient to give the following definitions: ∂U 1 ∂V , (5) ∇ζ · V h = + cos θ cos2 θ ∂λ ∂θ d ∂ 1 ∂ ∂ ∂ = + U + V cos θ + ζ˙ , (6) dt ∂t cos2 θ ∂λ ∂θ ∂ζ cos θ u cos θ v U= , V= , (7) a a 1 ∂ cos θ ∂ a ∇ζ = , , (8) cos θ a2 ∂λ a2 ∂θ where λ (longitude) and θ (latitude) are the coordinates. f is the Coriolis parameter, (u, v) are the wind components and (U, V) are the wind images. 2.2. Boundary conditions For each (Yin or Yang) subgrid LAM model, we impose Dirichlet boundary conditions in the horizontal and we define the top and bottom of the atmosphere to be material surfaces where the top is at constant pressure ptop and thus ζ˙ = 0 at

ζ = ζsurf , ζtop .

(9)

Given that the Yin and Yang grids are identical in size but not in rotation, it is therefore natural to give the same processor topology for both subdomains. In order to limit the amount of communication, the exchange of boundary conditions was optimized by allowing each individual processor for each subgrid to calculate in advance which points and processors of the other subgrid to send data to and/or receive data from. This is all tagged initially before the model starts to integrate. When there is a need to exchange boundary conditions, each processor interpolates and sends pre-assigned points to the receiving processor on the other subgrid. This limits the amount of data to be transmitted. OpenMP was also added to the interpolators to increase scalability. 2.3. Numerical solution ˆ e As in the global Canadian GEM operational model (Cot´ et al., 1998; Yeh et al., 2002), the same fully implicit semiLagrangian method is used in each subgrid of the Yin-Yang grid in order to solve dynamic core equations (1)–(4). This numerical method is based essentially on the following:

1915

right-hand sides of the above equations are then computed and added using the usual fractional step-time method (Yanenko, 1971). Consider a frictionless adiabatic prognostic equation of the form dF + G = 0, dt

(10)

where F represents one of the prognostic quantities Vh , ln (T/T ∗ ) − κBs, Bs + ln[1 + (∂B/∂ζ )s] and G represents the remaining terms, some of which are nonlinear. Such an equation is approximated by time differences and weighted averages along the trajectory determined by an approximate solution to dr = Vh (r, ζ , t), dt dζ = ζ˙ (r, ζ , t), dt

V2h d2 r = −r , dt 2 a2 d2 ζ = 0, dt 2

(11)

where r(λ, θ ) is the position vector on the local (to Yin or Yang) spherical coordinate system of the fluid element and a is the Earth’s radius. At each time step, the cubic Lagrange interpolation is used to update the static halo region of both panels of the Yin-Yang grid. Once the trajectory of the fluid element is known, that is at position (r, ζ ) and at forecast time t, then Eq. (10) may be integrated following the motion over a time interval t. Thus, F − F− 1 1 (12) + ( + )G + ( − )G− = 0. t 2 2 The scheme is decentred along the trajectory, as in Rivest et al. (1994). Cubic Lagrange interpolation is used everywhere for upstream evaluations (10) except for trajectory computations (11), where linear interpolation is used. In Eq. (12), grouping terms at the new time on the lefthand side and the known quantities on the right-hand side will result in the following system of nonlinear equations: F F− 1 − 2 − +G= −( )G , τ τ 1 + 2

(13)

where τ = (1 + 2)t/2. In this implicit treatment, the trajectory equation (11) is solved in a predictor–corrector manner (Yeh et al., 2002). 2.5. Spatial discretization

In each subdomain of the Yin-Yang grid, the vertical (1) two-time-level Crank–Nicholson iterative semi- discretization on a Charney–Phillips grid is as described Lagrangian method with time interpolation in an in Girard et al. (2010). Thus advecting wind; dVh ζ (2) iterative nonlinear solver for the Helmholtz problem (14) + f k × Vh + RT ∇ζ Bs + ∇ζ φ = 0, dt with a direct solver kernel; d T d ζ (3) iterative treatment of the Coriolis terms by grouping −κ ln B s + ζ˙ = 0, (15) them together with the nonlinear terms; dt T∗ dt (4) metric terms using the Lagrange multiplier approach d ζ Bs + ln 1 + δ B s ζ of Cˆot´e (1988). dt

2.4. Temporal discretization Equations (1)–(4) are first integrated in the absence of forcing, and the parametrized forcing terms appearing on the c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

ζ

+∇ζ · Vh + δζ ζ˙ + ζ˙ = 0, ∗ δ ζ − φ /RT ζ T = 0, − T∗ δζ (ζ + Bs) Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

(16) (17)

1916

A. Qaddouri and V. Lee

where the derivative is replaced by simple finite differences represented by the operator δζ and the averaging operators represented by overbars are introduced where it is required. The variables Vh , φ and B are on the same momentum vertical levels, whereas the variables T, ζ˙ are defined on thermodynamic levels that are staggered with respect to momentum levels. On each subgrid of the Yin-Yang grid, the horizontal spatial discretization is done as in Cˆot´e et al. (1998): a uniform discretization on a local Arakawa C grid is used. The discretization is centred and is then second-order in space. 2.6.

substructuring formulation of the elliptic problem and the use of a Krylov solver, as in Qaddouri (2008). The EBV solution starts with a vertical separation, where P is expanded in terms of vertical eigenvectors that ζ diagonalize the operator H = (δζ2 + δζ ). Equation (18) is then decoupled as the following set of Nk independent horizontal Helmholtz problems: ζ Pm + η(m)Pm = Rm H,

ζ P +

γ ζ (δ 2 P + δζ P ) = RE , κτ 2 RT ∗ ζ

η(m) =

2.7. Iterative formulation of the Schwarz method for an EBV problem In the optimized Schwarz method and as in the classical one, we solve the two discretized Helmholtz problems in the two subdomains l = 1, 2 of the Yin-Yang grid iteratively with iteration number k = 1, kmax : θj−1 cos2 1 (l),k Pi,j−1 + 2 P(l),k hλ cos2 θj i−1,j sin θj−1 sin θj−1 cos2 θj−1 2 1 − η+ 2 + ( 2 hλ cos θj sin θj−1 sin θj−1 2 cos θj 1 (l),k + ) Pi,j + 2 P(l),k sin θj hλ cos2 θj i+1,j cos2 θj + P(l),k = a2 Rh (l),k i,j , sin θj−1 sin θj i,j+1

(18)

θ˜j

RE is the right-hand side , Ctop and Csurf are constants and a 2 ζ =

1 ∂2 1 ∂ cos2 θ . + cos2 θ ∂λ2 ∂ sin θ ∂ sin θ

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

sin θ˜j sin θj

(19)

The above 3D EBV must be solved four times during each time step. We use here a domain decomposition method (DDM), where the solution of the global elliptic problem is obtained by iteratively solving the corresponding two subproblems separately on the Yin and the Yang grids and updating the values at the interfaces. This ‘classical’ DDM method corresponds to the method of Dirichlet interface conditions, which consists of using the subdomain solutions to update the interface conditions of neighbouring subdomains. To improve its performance we can include the solution derivatives and some precomputed optimized parameters in transmission conditions, as in Qaddouri et al. (2008). Further improvement can be made through a

(21)

i = 1, . . . , N; j = 1, . . . , M, where

ζ

(δζ Psurf + κPsurf ) = Csurf .

γ eig(m), m = 1, Nk, κτ 2 RT ∗

and eig(m) is the mth eingenvalue of the operator H. Because the DDM for the Yin-Yang grid is purely horizontal, we will concentrate on the DDM solution of the set of horizontal Helmholtz problems (20). Also, to make the equations more readable, the index m used to denote the number of the vertical mode will not appear in the following sections.

with the following top and bottom boundary conditions: δζ Ptop = Ctop ,

(20)

where Nk is the number of vertical levels, Rm H is the vertical projection of RE on the vertical mode number m,

Coupled nonlinear set of discretized equations

After spatial discretization, the resulting coupled set of nonlinear equations still has the form of Eq. (13). Terms on the right side, which involve upstream interpolation, are evaluated first. The coupled set is rewritten as a linear one (where the coefficients depend on the basic state) plus a nonlinear perturbation, which is placed on the right-hand side and which is relatively cheap to evaluate. The set is then solved iteratively by using the linear terms as a kernel and re-evaluating the nonlinear terms (including the Coriolis terms) on the right-hand side at each iteration using the most recent values of variables (horizontal wind components, . . . ). The nonlinearity is due mostly to logarithmic terms in the governing equations (1)–(4). The convergence and optimization of the iterative scheme were analyzed by Cˆot´e and Staniford (1988). Two nonlinear iterations, the minimum for stability, have been found sufficient for practical use. Since there are two outer iterations, the total number of iterations (and evaluation of the nonlinear terms) is four, giving a scheme that is very robust. The linear set of equations can be algebraically reduced to the solution of a 3D elliptic boundary value (EBV) problem for the composed variable P = φ + RT ∗ Bs, called generalized pressure, i.e.

m = 1, Nk,

θj + θj+1 , j = 1, M − 1, = 2 = sin θ˜j+1 − sin θ˜j , = sin θj+1 − sin θj ,

η is a constant eigenvalue and hλ = λi+1 − λi ; i = 1, N − 1 is the constant grid spacing along the longitude. At each iteration k we use the following discretizations of the Dirichlet transmission conditions on each interface

d(l) (d = 1, · · · , 4): P(l),k = P(3−l),k−1 . 3.

(22)

Physical parametrization

The forcing terms appearing in the right-hand sides of Eqs (1) and (2) are specified or parametrized by the unified Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

(a)

(b)

(c)

(d)

1917

(e)

Figure 2. Objective evaluation of 120 h forecasts against radiosondes for the GEM global (blue lines) and Yin-Yang (red lines) global forecast systems for (a) U winds in m s−1 , (b) wind speed in m s−1 , (c) geopotential height (dam), (d) temperature (K) and (e) dew point depression (K). The RMSE (solid lines) and bias (dashed lines) are shown. Verification is performed over North America, for a set of 42 integrations over the winter period 15 December 2008–14 February 2009.

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1918

A. Qaddouri and V. Lee

(a)

(b)

(c)

(d)

(e)

Figure 3. Objective evaluation of 120 h forecasts against radiosondes for the GEM global (blue lines) and Yin-Yang (red lines) global forecast systems for (a) U winds in m s−1 , (b) wind speed in m s−1 , (c) geopotential height (dam), (d) temperature (K) and (e) dew point depression (K). The RMSE (solid lines) and bias (dashed lines) are shown. Verification is performed over North America using KF trigger=0.05, for a set of 42 integrations over the summer period 18 June 2008–18 August 2008.

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

(a)

(b)

(c)

(d)

1919

(e)

Figure 4. Objective evaluation of 120 h forecasts against radiosondes for the GEM global (blue lines) and Yin-Yang (red lines) global forecast systems for (a) U winds in m s−1 , (b) wind speed in m s−1 , (c) geopotential height (dam), (d) temperature (K) and (e) dew point depression (K). The RMSE (solid lines) and bias (dashed lines) are shown. Verification is performed over North America using KF trigger=0.0475, for a set of 42 integrations over the summer period 18 June 2008–18 August 2008.

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1920

A. Qaddouri and V. Lee

(a)

(b)

(c)

(d)

Figure 5. 24 h accumulated precipitation for integration using analysis on 11 August 2010 at 0000 UTC. (a) GEM global (uniform grid), (b) Yin-Yang, (c) Yin-Yang (KF trigger=0.0475) and (d) GEM global (variable grid) with Yin grid in core, all using KF trigger=0.05 except for (c).

(9) the ozone climatology based on ozonesonde and RPN (Recherche Prevision Numerique) physics package. The choice of parametrization depends on the space- and satellite measurements from Fortuin and Kelder time-scales of the forecast. In this project, the physical (Fortuin and Kelder, 1998). processes are the same as in the GEM global forecast system even though we are coupling two LAM grids (the two panels 4. Digital filter and diffusion of the Yin-Yang grid). The parametrization scheme used for the physics configuration is as follows: For both the GEM global and Yin-Yang runs, the diabatic (1) the Interaction Soil Biosphere Atmosphere (ISBA) digital filter (Fillion et al., 1995) with a 6 hour span land-surface scheme for the surface-layer effects was applied at the beginning of the model integration to dampen high-frequency noise from the initial conditions. (B´elair et al., 2003a,b); (2) Geleyn boundary-layer cloud scheme (Geleyn, 1987); In both models, horizontal diffusion is expressed in the (3) Kuo transient shallow convection scheme (Kuo, form of a scale-selective hyper-Laplacian (∇ 6 ) applied to 1974); momentum variables and temperature. This hyper-diffusion (4) the Kain–Fritsch deep convection scheme (Kain and is solved by an implicit scheme described in Qaddouri and Fritsch, 1993); Lee (2008). (5) the Bougeault–Lacarr`ere turbulent mixing-length An enhanced diffusion is applied in the sponge layer scheme (Bougeault and Lacarr`ere, 1989); covering the six uppermost levels. In the GEM global run, (6) the radiative transfer scheme from Li and Barker (Li it uses an implicit horizontal (∇ 2 ) scheme, whereas for the and Barker, 2005); (7) the non-orographic gravity-wave drag scheme by Yin-Yang runs only the explicit nine-point filter scheme was available under the GEM–LAM mode. The coefficient in the Hines (Hines, 1997a,b); (8) the inclusion of a methane oxidation parametrization latter scheme was adjusted to give the closest behaviour to scheme (same scheme used in the ECMWF model); the implicit scheme. c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

1921

(a)

(b)

Figure 6. Quantitative precipitation forecasts for 24–48 h for the GEM global (blue lines with crosses) and Yin-Yang (red lines with stars) global forecasts. The bias is calculated in such a way that when the value is closer to 1.0 there is less bias. The higher threat score indicates more precision in the forecasts. Verification is carried out over North America using the observational surface network from the Standard Hydrometeorological Exchange Format(SHEF) for a set of 42 integrations over the winter period 15 December 2008–14 February 2009.

5. Experimental set-up and evaluation

5.1. Objective evaluation

The goal in this preliminary objective evaluation is to demonstrate that we could use a Yin-Yang grid in the GEM model to produce the same forecast results as the GEM global configuration. The latter uses a global latitude–longitude grid. In order to make a fair comparison in terms of having the closest horizontal resolution, the GEM global had to be configured with twice the number of points along the longitude versus the number of points along the latitude. Our current 33 km operational GEM global has a horizontal grid of 800 × 600 points with 80 vertical levels. We decided to compare a global set-up using a horizontal grid of 800 × 400 points and 80 vertical levels with the Yin-Yang grid system of (600 × 200) × 2 points and 80 vertical levels. Geophysical fields were generated for each respective grid: the GEM global using 800 × 400 grid points and 600 × 200 grid points on each Yin and Yang grid.

For each system (Yin-Yang and GEM global), we conducted two series of tests: 42 winter cases and 42 summer cases where each case is a five-day forecast verified against upper-air and radiosonde observations. Each season spans over a two-month period at every 36 hours. The values of the quantitative precipitation forecasts (QPFs) were also compared between the two systems. This is a typical method used to evaluate models at CMC. When we first compared the objective evaluation of both systems for the winter series against upper-air and radiosonde observations, both gave very similar results, as shown in Figure 2. The blue lines represent the GEM global runs, which are our control runs, and the red lines represent the Yin-Yang runs. Note that the bias and the root-meansquare errors (RMSE) are almost equivalent for both the GEM global runs and the Yin-Yang runs. Each result from each system is almost superimposed for the winds (UU, UV),

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1922

A. Qaddouri and V. Lee

(a)

(b)

Figure 7. Quantitative precipitation forecasts for 24–48 h for the GEM global (blue lines with crosses) and Yin-Yang (red lines with stars) global forecast. The bias is calculated in such a way that when the value is closer to 1.0 there is less bias. The higher threat score indicates more precision in the forecasts. Verification is carried out over North America using the observational surface network from the Standard Hydrometeorological Exchange Format (SHEF) for a set of 42 integrations over the summer period 18 June 2008–18 August 2008.

temperature (TT), geopotential height (GZ) and dewpoint depression (ES). We then ran the summer series, where the results in Figure 3 showed a slight unfavourable significant difference for the Yin-Yang system from the bias in the GZ above 100 mb. On a closer inspection, we noticed that the YinYang GEM model was underpredicting larger amounts of precipitation, particularly for rainfall over 20 mm in a 24 h period and specifically over land. It was also noted that the vertical motion was a bit less intense (smoother) in the Yin-Yang than in the GEM global configuration. It is known that the intensity of the vertical motion is directly related to the resolution of the grid. The factor in the Kain–Fritsch (KF) deep convection scheme is usually adjusted to suit the average resolution of the grid in order to control how easily the KF scheme can be triggered. The finer the resolution, the higher the grid-scale vertical velocity. When the temperature of a buoyant parcel plus temperature perturbation associated with the grid-scale vertical velocity c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

is higher than the environmental temperature, then the KF scheme is triggered. Unfortunately, in our current operational models the KF trigger value is tuned to be dependent on an averaged horizontal resolution and not on the local resolution nor on the true latitude. We found that lowering the KF trigger value slightly to 0.0475 m s−1 was sufficient to obtain the same validation as the GEM global for the summer series (Figure 4). Note that there is no longer a significant difference in the bias of the GZ, as well as there being a slight improvement in the bias of TT (temperature) between 400 and 250 mb. The winter series for the Yin-Yang model were repeated with the newly adjusted KF trigger value of 0.0475 m s−1 and this had no impact on the results. One of the biggest challenges in comparing the two systems is to match the grids as closely as possible. The Yin-Yang grid remains quasi-uniform while the global grid has an increasing resolution for grid points approaching the poles. The grids in GEM are designed in such a way that Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

1923

(a)

(b)

Figure 8. The 72 h forecast of the air temperature (TT) close to the surface for (a) GEM global (uniform grid) and (b) Yin-Yang grid.

even with almost the same resolution at the Equator, the GEM global grid points are staggered to the grid points of the Yin-Yang grid. To test the sensitivity of precipitation in the location of grid points, we decided to simulate a global run using the Yin grid as the core domain in a global variable grid described in Cˆot´e et al. (1993). Figure 5 shows four plots of the same region of the USA with a 24 h forecast snapshot of accumulated precipitation. All use a KF trigger of 0.05 m s−1 except for Figure 5(c). Note that the precipitation maxima of the global uniform grid shown in Figure 5(a) are much greater than those in the original Yin-Yang run shown in Figure 5(b). In Figure 5(c), the maxima of the precipitation increased for the Yin-Yang grid, due to decreasing the KF trigger to 0.0475 m s−1 . In Figure 5(d), the precipitation shown is from a run using the global variable grid using the Yin grid as its core uniform domain and it also used the original KF trigger of 0.05 m s−1 . Note that Figure 5(d) bore the greatest resemblance to the original Yin-Yang run in Figure 5(b) because the grid points were the most closely matched. This indicates that the precipitation pattern from the GEM can be affected not only by the horizontal resolution of grid points but also by their location. One of the ways of validating precipitation at CMC is to look at the quantitative precipitation forecasts (QPFs) from each model of both the winter cases (Figure 6) and the summer cases (Figure 7). In each figure there are two c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

plots: one is the bias and the other is the threat score, where the blue represents GEM global and the red represents YinYang. The bias calculated is based on the number of forecasts for a category divided by the number of observations in this category. The threat score is a measure of the number of correct forecasts for a category divided by the number of events forecast or observed in this category. In looking at these two QPF plots, it seems that there are no significant meteorological differences. A major advantage in using the Yin-Yang grid system to represent the Earth’s spherical geometry is that one avoids singularities and mesh convergence in the polar regions. There is an expectation that a Yin-Yang model may give better results in the regions near the poles. In all the integrations done with the GEM global model, we can find numerical noise in the pole areas and an example is shown in Figure 8(a) for the 72 h forecast of air temperature (TT) close to the surface. The same forecast produced by the YinYang model shows much smoother contours at the south pole and maintains the same pattern as the GEM global elsewhere (see Figure 8(b)). 5.2. Preliminary parallel performance As in the experimental set-up, we have set the resolution to be the same at the Equator for both systems (Yin-Yang and GEM global). We have chosen to use fast Fourier transforms Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1924

A. Qaddouri and V. Lee

Figure 9. Wall-clock time in seconds and speed-up versus number of CPUs for 96 time steps. Both GEM global and Yin-Yang runs are configured for 39 km resolution at the Equator.

(FFT) in the EBV solver for both systems, which meant that the choice of the number of longitude points (N) is limited. In the case of Yin-Yang, there is an additional constraint whereby the number of longitude points (N) must be three times the number of latitude points (M). For the global grid, N is twice M. The resolution chosen for this preliminary parallel performance test is 39 km where the global grid is 1024 × 512 and each Yin-Yang subgrid is 799 × 267 points. Both have the same definition for the 80 vertical levels and both use the same configuration for the physics. The time step given is 720 s. For each test, we let the model run 96 time steps and measured the wall-clock time. We tried various processor topologies for the two systems and retained the best ones for each system. Our best approach to compare performance, especially in the context of operational runs, is to plot the wall-clock time versus the number of CPUs used in each run. In Figure 9(a), we see plots of both systems with the wall-clock time of a run versus the number of CPUs used in a given run. The red plot line represents the Yin-Yang runs, and this starts to outperform the GEM global run (blue line) when using more than 512 CPUs. Another way to compare the parallel scalability of the two models is to plot the parallel speed-up versus the number of CPUs used. In Figure 9(b), the speed-up of 1.0 is based on the runs c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

using 192 CPUs, as we cannot run a 39 km resolution model on 1 CPU. The results plotted are the best timings for the number of CPUs given to each system, since it is rare that both can use the same number of CPUs. This figure shows a better overall scalability with the Yin-Yang model. Looking at Table 1, we can see that the advection scales better in the Yin-Yang case because of the LAM mode when we are increasing the number of processors. At the same time, the smaller value in the physics computation indicates that there are fewer points to compute for the physical processes. The disadvantage of Yin-Yang is in the computation for the solver, as it must solve iteratively between Yin and Yang to obtain its solution. The other additional cost in the YinYang system is the extra communication between the two panels in the update of the boundary conditions and other intermediate steps to complete the two-way coupling. The code still needs further optimization and we hope to see an even better speed-up in the future. 6.

Conclusion

The future global forecasting Yin-Yang model has been developed by a two-way coupling method between two GEM LAMs discretized on the two panels of the Yin-Yang grid. A Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

1925

Table I. Performance breakdown of timings in seconds. Mode Yin-Yang Global Yin-Yang Global Yin-Yang Global Yin-Yang Global

PE Topo

NEST

BAC

SOL

ADV

PHY

TOT

CPUs

(12 × 8 × 1) × 2 (3 × 32 × 2) (20 × 8 × 1) × 2 (2 × 40 × 4) (16 × 12 × 1) × 2 (3 × 32 × 4) (16 × 8 × 2) × 2 (4 × 37 × 4)

64 n/a 54 n/a 43 n/a 32 n/a

38 6 19 5 19 4 16 4

122 39 125 35 73 28 67 23

134 162 93 127 74 109 65 82

233 272 130 222 113 145 90 115

676 597 496 508 386 388 326 325

192 192 320 320 384 384 512 592

PE Topo –processor topology (Npex × Npey × OpenMP). NEST –exchange of boundary conditions. BAC –Back Substitution. SOL –Solver. ADV –Advection. PHY –Physics.TOT –Total Wallclock. CPUs –Total Number of CPUs used.

comparison has been made by matching as closely as possible the grid resolution of each system (Yin-Yang and GEM global) and then performing a model validation for both systems. Preliminary results of the new Yin-Yang system show that it is capable of producing meteorological forecasts equivalent to the GEM global model. One big advantage in using this type of system is to avoid the singularities and convergence of meridians in the polar regions. Another benefit is to have a more balanced computational load for scalability purposes. Our next step would be to compare objective evaluations with a much higher resolution grid and using a non-hydrostatic version of the GEM model. The results for parallel performance also indicate better scalability compared with the global model when the system is given more processors. Investigation and development will continue, to consider whether this Yin-Yang system could become the future global forecasting system at CMC. Acknowledgement ˆ e and Claude The authors thank their colleagues Jean Cot´ Girard and the two anonymous reviewers for their helpful comments. References Adcroft A, Campin J-M, Hill C, Marshal J. 2004. Implementation of atmosphere–ocean general circulation model on the spherical cube. Mon. Weather Rev. 132: 2845–2863. Baba Y, Takahashi K, Sugimura T. 2010. Dynamical core of an atmospheric general circulation model, on Yin-Yang grid. Mon. Weather Rev. 138(10): 3988–4005, DOI:10.1175/2010MWR3375.1. B´elair S, Crevier L-P, Mailhot J, Bilodeau B, Delage Y. 2003a. Operational implementation of the ISBA land surface scheme in the Canadian regional weather forecast model. Part I: Warm season results. J. Hydrometeorol. 4: 352–370. B´elair S, Brown R, Mailhot J, Bilodeau B, Crevier L-P. 2003b. Operational implementation of the ISBA land surface scheme in the Canadian regional weather forecast model. Part II: Cold season results. J. Hydrometeorol. 4: 371–386. B´elair S, Roch M, Leduc A-M, Vaillancourt P-A, Laroche S, Mailhot J. 2009. Medium-range quantitative precipitation forecasts from Canada’s new 33 km deterministic global operational system. Weather and Forecasting 24: 690–708. Bougeault P, Lacarr`ere P. 1989. Parameterization of orography-induced turbulence in a mesobeta-scale model. Mon. Weather Rev. 117: 1872–1890. Cˆot´e J. 1988. A Lagrange multiplier approach for the metric terms of semi-Lagrangian models on the sphere. Q. J. R. Meteorol. Soc. 114: 1347–1352. Cˆot´e J, Staniforth A. 1988. A two-time-level semi-Lagrangian semiimplicit scheme for spectral models. Mon. Weather Rev. 116: 2003–2012. c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Cˆot´e J, Roch M, Staniforth A, Fillion L. 1993. A variable-resolution semi-Lagrangian finite-element global model of the shallow-water equations. Mon. Weather Rev. 121: 231–243. Cˆot´e J, Gravel S, M´ethot A, Patoine A, Roch M, Staniforth A. 1998. The operational CMC–MRB Global Environmental Multiscale (GEM) model. Part I: Design considerations and formulation. Mon. Weather Rev. 126: 1373–1395. Dudhia J, Bresch J-F. 2002. A global version of PSU-NCAR mesoscale model. Mon. Weather Rev. 130: 2989–3007. Fillion L, Mitchel H-L, Ritchie H-R, Staniforth A. 1995. The impact of a digital filter finalization technique in a global data assimilation system. Tellus 47A: 304–323. Fortuin P, Kelder H. 1998. An ozone climatology based on ozonesonde and satellite measurements. J. Geophys. Res. 103: 31709–31734. Geleyn J-F. 1987. ‘Use of a modified Richardson number for parameterizing the effect of shallow convection. Short- and medium-range numerical weather prediction’. In WMO/IUGG NWP Symposium. WMO: Tokyo, Japan; pp 141–149. Girard C, Plante A, Gravel S, Qaddouri A, Chamberland A, Spacek L, Lee A, Desgagn´e M. 2010. GEM4.1: A non-hydrostatic atmospheric model (Euler equations). RPN document. http://collaboration. cmc.ec.gc.ca/science/rpn/publications/pdf/GEM4.1.pdf. Hines C-O. 1997a. Doppler-spread parameterization of gravity-wave momentum deposition in the middle atmosphere. Part 1: Basic formulation. J. Atmos. Solar–Terr. Phys. 59: 371–386. Hines C-O. 1997b. Doppler-spread parameterization of gravity-wave momentum deposition in the middle atmosphere. Part 2: Broad and quasi monochromatic spectra, and implementation. J. Atmos. Solar–Terr. Phys. 59: 387–400. Kageyama A, Sato T. 2004. The ‘Yin-Yang grid’: An overset grid in spherical geometry. Geochem. Geophys. Geosystems 5: Q09005. DOI:10.1029/2004GC000734. Kain J, Fritsch J-M. 1993. Convective parameterization for mesoscale models: The Kain–Fritsch scheme. The representation of cumulus convection in numerical models. Meteorol. Monogr. 46: 165–170. Kuo H-L. 1974. Further studies of the parameterization of the influence of cumulus convection on large-scale flow. J. Atmos. Sci. 31: 1232–1240. Li J, Barker H-W. 2005. A radiation algorithm with correlated-k distribution. Part I: Local thermal equilibrium. J. Atmos. Sci. 62: 286–309. Qaddouri A. 2008. Optimized Schwarz methods with the Ying-Yang grid for shallow water equations. In Domain Decomposition Methods in Science and Engineering XVII. Lecture Notes in Computational Science and Engineering 60 347–353. Qaddouri A. 2011. Nonlinear shallow-water equations on the Yin-Yang grid. Q. J. R. Meteorol. Soc. 137: 810–818, DOI:10.1002/qj.792. Qaddouri A, Lee V. 2008. ‘Solution of the implicit formulation of high order diffusion for the Canadian Atmospheric GEM model’. In Proc. 2008 Spring Simulation Multiconf., High Performance Computing Symp., Hamilton JA Jr et al. (eds). Soc. For Modeling and Simulation Internat.: Ottawa, Canada; pp 362–367. Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1926

A. Qaddouri and V. Lee

Qaddouri A, Laayouni L, Loisel S, Cˆot´e J, Gander M-J. 2008. Optimized Schwarz methods with an overset grid for the shallowwater equations: Preliminary results. Appl. Numerical Math. 58: 459–471. Rivest C, Staniforth A, Robert A. 1994. Spurious resonant response of semi-Lagrangian discretizations to orographic forcing: Diagnosis and solution. Mon. Weather Rev. 122: 366–376.

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Yanenko N-N. 1971. The method of Fractional Steps. Springer: Berlin; 160 pp. Yeh K-S, Cˆot´e J, Gravel S, M´ethot A, Patoine A, Roch M, Staniforth A. 2002. The CMC–MRB global environmental multiscale (GEM) model: Part III – Nonhydrostatic formulation. Mon. Weather Rev. 130: 339–356.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Q. J. R. Meteorol. Soc. 137: 1913–1926, October 2011 A

The Canadian Global Environmental Multiscale model on the Yin-Yang grid system Abdessamad Qaddouri* and Vivian Lee Meteorological Research Division, Atmospheric Science and Technology Directorate, Environment Canada, Dorval, Qu´ebec, Canada *Correspondence to: A. Qaddouri, Environment Canada, 2121 Trans-Canada Hwy, Dorval QC H9P 1J3, Canada. E-mail: [email protected]

At the Canadian Meteorological Center (CMC), we are currently developing the future global forecasting Yin-Yang model. In the horizontal we use spherical coordinates on the overset Yin-Yang grid, while in the vertical we use a loghydrostatic-pressure coordinate on the Charney–Phillips grid. The parametrization of physical processes is kept the same as in the current Global Environmental Multiscale (GEM) operational model. The Yin-Yang global forecast is performed by considering a domain decomposition (a two-way coupling method) between two limited-area models (LAMs) discretized on the two panels of the Yin-Yang grid and using the same time step. Each panel of the Yin-Yang grid system is extended by a static halo region and uses the same fully implicit semi-Lagrangian method as in the GEM operational model to solve its own dynamic core. The spatial and time discretizations are implemented independently on each quasi-uniform latitude–longitude subgrid. The static halo region plays the same role as the piloting region in limited-area modelling. Since the two subgrids of the Yin-Yang grid do not match, the update of the variables in the pilot region is done by cubic Lagrange interpolation. For our model validation, we ran 42 winter and 42 summer cases using analysis from 2008–2009 and we compared five-day forecast results against observations. No noise is seen in the overlap regions during the simulations. Preliminary results presented in this article are encouraging and demonstrate that in comparison with observations the new Yin-Yang system performs as well as the GEM global c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd. model. Key Words: method

overset grid on the sphere; semi-Lagrangian semi-implicit scheme; domain decomposition

Received 22 December 2010; Revised 3 May 2011; Accepted 25 May 2011; Published online in Wiley Online Library 1 August 2011 Citation: Qaddouri A, Lee V. 2011. The Canadian Global Environmental Multiscale model on the Yin-Yang grid system. Q. J. R. Meteorol. Soc. 137: 1913–1926. DOI:10.1002/qj.873

1. Introduction The deterministic Global Environmental Multiscale (GEM) operational model (Yeh et al., 2002) constitutes a very important tool for medium-range numerical weather prediction at the Canadian Meteorological Center (CMC). Over the last decade, significant improvements have been made in the GEM physical process representation (B´elair

et al., 2009) and the model lid has been raised to 0.1 mb with no major changes in its dynamical core part. In the current year, however, the non-staggered vertical grid with a hydrostatic pressure coordinate in the GEM model (Yeh et al., 2002) was replaced by the Charney–Phillips grid with a log-hydrostatic-pressure coordinate (Girard et al., 2010). This new spatial vertical discretization has a positive impact on model performance, since it improves the accuracy of

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

1914

A. Qaddouri and V. Lee

area modelling (LAM). Each subdomain is extended by a static halo region that plays the same role as the piloting region in the LAM. The Schwarz method consists of solving the subdomain problems iteratively by using the subdomain solutions to update the boundary conditions in the pilot region of the neighbouring subdomain. Since the two subgrids of the Yin-Yang grid do not match, a full cubic Lagrange interpolation is used to bring the boundary conditions to each respective target grid. This technique has been used to solve the nonlinear shallow-water equations on the Yin-Yang grid in Qaddouri (2011). The three-dimensional (3D) dynamical core has already been developed and tested on both the polar-cap grid (Dudhia and Bresch, 2002) and the cubed sphere grid (Adcroft et al., 2004). The Yin-Yang grid has just been considered for the first time in the development of a 3D dynamical core for simulating analytical cases of atmospheric flow ( Baba et al., 2010). In this article, we will focus on the feasibility of the Yin-Yang grid for operational medium-range global forecasts while considering both dynamical and physical process representations. The model formulation and numerical solution are described in section 2. For simplicity, only hydrostatic model formulations will be shown, as the non-hydrostatic equations involve a few extra terms but do not change the technique used in obtaining the numerical solution. The schemes used for the parametrization of physical processes are presented in section 3. The use of a digital filter and horizontal diffusion is described in section 4. The numerical results are shown and discussed in section 5. Finally, our conclusions are presented in section 6.

(a)

(b)

Figure 1. (a) Latitude–longitude grid and (b) Yin-Yang grid. This figure is available in colour online at wileyonlinelibrary.com/journal/qj

2. 2.1.

the hydrostatic and non-hydrostatic equations by defining the vertical finite differences at logarithmic midpoints. This staggered vertical grid GEM version is expected to become operational in 2011, and we will refer to it as the GEM global model. A problem associated with the construction of any global atmospheric model is the numerical representation of the spherical geometry. Currently, the GEM global model uses the latitude–longitude (denoted Lat–Lon) grid system, which leads to singularities and convergence of meridians in the polar regions. Using today’s computer architectures with distributed memory, this may even result in an unbalanced computational load in the context of domain decomposition with message-passing interfaces. In this paper, we developed a Yin-Yang model by replacing a Lat–Lon mesh grid as in Figure 1(a) with a Yin-Yang grid system (Figure 1(b)) that imposes a more uniform grid spacing. The Yin-Yang model uses the Schwarz domain decomposition method (DDM) in order to solve a system of forced primitive equations (PEs) on a quasi-uniform overset grid free from polar singularity. This grid system is constructed by overlapping two perpendicularly oriented identical parts of a global Lat–Lon grid (Figure 1(b)) as suggested by Kageyama and Sato (2004). In our implementation of the DDM, the governing equations on each (Yin or Yang) subgrid are the PEs and these are discretized using the implicit semi-Lagrangian method as in the GEM operational model (Cˆot´e et al., 1998; Yeh et al., 2002). The solution for each subdomain is based on the same technique used in limited c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Model formulation Governing equations

In each of the two subdomains of the Yin-Yang grid system, the governing equations are the forced hydrostatic primitive equations described in Girard et al. (2010) as dVh + f k × Vh + RT∇ζ Bs + ∇ζ φ = FH , dt T d d ˙ ln −κ (Bs) + ζ = F T , dt T∗ dt ∂B d ∂ ζ˙ Bs + ln 1 + s + ∇ζ · Vh + + ζ˙ = 0, dt ∂ζ ∂ζ ∗ ∂ ζ − φ /RT T = 0, − T∗ ∂ (ζ + Bs)

(1) (2) (3)

(4)

where ζ = ln p − Bs is a log-hydrostatic-pressure type coordinate with p the pressure, B = [(ζ − ζtop )/(ζsurf − ζtop )] the metric term, s = ln psurf − ζsurf the mass variable that depends only on the horizontal, T the temperature, φ = φ ∗ + φ the geopotential, φ ∗ (ζ ) = −RT ∗ (ζ − ζsurf ), ∗ T = const and Vh the horizontal wind. R is the gas constant and κ = R/cp , where cp is the specific heat at constant pressure. The coordinate system permits a switch-controlled choice between the hydrostatic primitive equations and non-hydrostatic Euler equations. See Girard et al. (2010) for additional details. The prognostic equations (1)–(3) are respectively the horizontal momentum, thermodynamic and continuity Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

equations, and (4) is the diagnostic hydrostatic equation. The terms FH and F T are parametrized physical forcings. It is also convenient to give the following definitions: ∂U 1 ∂V , (5) ∇ζ · V h = + cos θ cos2 θ ∂λ ∂θ d ∂ 1 ∂ ∂ ∂ = + U + V cos θ + ζ˙ , (6) dt ∂t cos2 θ ∂λ ∂θ ∂ζ cos θ u cos θ v U= , V= , (7) a a 1 ∂ cos θ ∂ a ∇ζ = , , (8) cos θ a2 ∂λ a2 ∂θ where λ (longitude) and θ (latitude) are the coordinates. f is the Coriolis parameter, (u, v) are the wind components and (U, V) are the wind images. 2.2. Boundary conditions For each (Yin or Yang) subgrid LAM model, we impose Dirichlet boundary conditions in the horizontal and we define the top and bottom of the atmosphere to be material surfaces where the top is at constant pressure ptop and thus ζ˙ = 0 at

ζ = ζsurf , ζtop .

(9)

Given that the Yin and Yang grids are identical in size but not in rotation, it is therefore natural to give the same processor topology for both subdomains. In order to limit the amount of communication, the exchange of boundary conditions was optimized by allowing each individual processor for each subgrid to calculate in advance which points and processors of the other subgrid to send data to and/or receive data from. This is all tagged initially before the model starts to integrate. When there is a need to exchange boundary conditions, each processor interpolates and sends pre-assigned points to the receiving processor on the other subgrid. This limits the amount of data to be transmitted. OpenMP was also added to the interpolators to increase scalability. 2.3. Numerical solution ˆ e As in the global Canadian GEM operational model (Cot´ et al., 1998; Yeh et al., 2002), the same fully implicit semiLagrangian method is used in each subgrid of the Yin-Yang grid in order to solve dynamic core equations (1)–(4). This numerical method is based essentially on the following:

1915

right-hand sides of the above equations are then computed and added using the usual fractional step-time method (Yanenko, 1971). Consider a frictionless adiabatic prognostic equation of the form dF + G = 0, dt

(10)

where F represents one of the prognostic quantities Vh , ln (T/T ∗ ) − κBs, Bs + ln[1 + (∂B/∂ζ )s] and G represents the remaining terms, some of which are nonlinear. Such an equation is approximated by time differences and weighted averages along the trajectory determined by an approximate solution to dr = Vh (r, ζ , t), dt dζ = ζ˙ (r, ζ , t), dt

V2h d2 r = −r , dt 2 a2 d2 ζ = 0, dt 2

(11)

where r(λ, θ ) is the position vector on the local (to Yin or Yang) spherical coordinate system of the fluid element and a is the Earth’s radius. At each time step, the cubic Lagrange interpolation is used to update the static halo region of both panels of the Yin-Yang grid. Once the trajectory of the fluid element is known, that is at position (r, ζ ) and at forecast time t, then Eq. (10) may be integrated following the motion over a time interval t. Thus, F − F− 1 1 (12) + ( + )G + ( − )G− = 0. t 2 2 The scheme is decentred along the trajectory, as in Rivest et al. (1994). Cubic Lagrange interpolation is used everywhere for upstream evaluations (10) except for trajectory computations (11), where linear interpolation is used. In Eq. (12), grouping terms at the new time on the lefthand side and the known quantities on the right-hand side will result in the following system of nonlinear equations: F F− 1 − 2 − +G= −( )G , τ τ 1 + 2

(13)

where τ = (1 + 2)t/2. In this implicit treatment, the trajectory equation (11) is solved in a predictor–corrector manner (Yeh et al., 2002). 2.5. Spatial discretization

In each subdomain of the Yin-Yang grid, the vertical (1) two-time-level Crank–Nicholson iterative semi- discretization on a Charney–Phillips grid is as described Lagrangian method with time interpolation in an in Girard et al. (2010). Thus advecting wind; dVh ζ (2) iterative nonlinear solver for the Helmholtz problem (14) + f k × Vh + RT ∇ζ Bs + ∇ζ φ = 0, dt with a direct solver kernel; d T d ζ (3) iterative treatment of the Coriolis terms by grouping −κ ln B s + ζ˙ = 0, (15) them together with the nonlinear terms; dt T∗ dt (4) metric terms using the Lagrange multiplier approach d ζ Bs + ln 1 + δ B s ζ of Cˆot´e (1988). dt

2.4. Temporal discretization Equations (1)–(4) are first integrated in the absence of forcing, and the parametrized forcing terms appearing on the c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

ζ

+∇ζ · Vh + δζ ζ˙ + ζ˙ = 0, ∗ δ ζ − φ /RT ζ T = 0, − T∗ δζ (ζ + Bs) Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

(16) (17)

1916

A. Qaddouri and V. Lee

where the derivative is replaced by simple finite differences represented by the operator δζ and the averaging operators represented by overbars are introduced where it is required. The variables Vh , φ and B are on the same momentum vertical levels, whereas the variables T, ζ˙ are defined on thermodynamic levels that are staggered with respect to momentum levels. On each subgrid of the Yin-Yang grid, the horizontal spatial discretization is done as in Cˆot´e et al. (1998): a uniform discretization on a local Arakawa C grid is used. The discretization is centred and is then second-order in space. 2.6.

substructuring formulation of the elliptic problem and the use of a Krylov solver, as in Qaddouri (2008). The EBV solution starts with a vertical separation, where P is expanded in terms of vertical eigenvectors that ζ diagonalize the operator H = (δζ2 + δζ ). Equation (18) is then decoupled as the following set of Nk independent horizontal Helmholtz problems: ζ Pm + η(m)Pm = Rm H,

ζ P +

γ ζ (δ 2 P + δζ P ) = RE , κτ 2 RT ∗ ζ

η(m) =

2.7. Iterative formulation of the Schwarz method for an EBV problem In the optimized Schwarz method and as in the classical one, we solve the two discretized Helmholtz problems in the two subdomains l = 1, 2 of the Yin-Yang grid iteratively with iteration number k = 1, kmax : θj−1 cos2 1 (l),k Pi,j−1 + 2 P(l),k hλ cos2 θj i−1,j sin θj−1 sin θj−1 cos2 θj−1 2 1 − η+ 2 + ( 2 hλ cos θj sin θj−1 sin θj−1 2 cos θj 1 (l),k + ) Pi,j + 2 P(l),k sin θj hλ cos2 θj i+1,j cos2 θj + P(l),k = a2 Rh (l),k i,j , sin θj−1 sin θj i,j+1

(18)

θ˜j

RE is the right-hand side , Ctop and Csurf are constants and a 2 ζ =

1 ∂2 1 ∂ cos2 θ . + cos2 θ ∂λ2 ∂ sin θ ∂ sin θ

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

sin θ˜j sin θj

(19)

The above 3D EBV must be solved four times during each time step. We use here a domain decomposition method (DDM), where the solution of the global elliptic problem is obtained by iteratively solving the corresponding two subproblems separately on the Yin and the Yang grids and updating the values at the interfaces. This ‘classical’ DDM method corresponds to the method of Dirichlet interface conditions, which consists of using the subdomain solutions to update the interface conditions of neighbouring subdomains. To improve its performance we can include the solution derivatives and some precomputed optimized parameters in transmission conditions, as in Qaddouri et al. (2008). Further improvement can be made through a

(21)

i = 1, . . . , N; j = 1, . . . , M, where

ζ

(δζ Psurf + κPsurf ) = Csurf .

γ eig(m), m = 1, Nk, κτ 2 RT ∗

and eig(m) is the mth eingenvalue of the operator H. Because the DDM for the Yin-Yang grid is purely horizontal, we will concentrate on the DDM solution of the set of horizontal Helmholtz problems (20). Also, to make the equations more readable, the index m used to denote the number of the vertical mode will not appear in the following sections.

with the following top and bottom boundary conditions: δζ Ptop = Ctop ,

(20)

where Nk is the number of vertical levels, Rm H is the vertical projection of RE on the vertical mode number m,

Coupled nonlinear set of discretized equations

After spatial discretization, the resulting coupled set of nonlinear equations still has the form of Eq. (13). Terms on the right side, which involve upstream interpolation, are evaluated first. The coupled set is rewritten as a linear one (where the coefficients depend on the basic state) plus a nonlinear perturbation, which is placed on the right-hand side and which is relatively cheap to evaluate. The set is then solved iteratively by using the linear terms as a kernel and re-evaluating the nonlinear terms (including the Coriolis terms) on the right-hand side at each iteration using the most recent values of variables (horizontal wind components, . . . ). The nonlinearity is due mostly to logarithmic terms in the governing equations (1)–(4). The convergence and optimization of the iterative scheme were analyzed by Cˆot´e and Staniford (1988). Two nonlinear iterations, the minimum for stability, have been found sufficient for practical use. Since there are two outer iterations, the total number of iterations (and evaluation of the nonlinear terms) is four, giving a scheme that is very robust. The linear set of equations can be algebraically reduced to the solution of a 3D elliptic boundary value (EBV) problem for the composed variable P = φ + RT ∗ Bs, called generalized pressure, i.e.

m = 1, Nk,

θj + θj+1 , j = 1, M − 1, = 2 = sin θ˜j+1 − sin θ˜j , = sin θj+1 − sin θj ,

η is a constant eigenvalue and hλ = λi+1 − λi ; i = 1, N − 1 is the constant grid spacing along the longitude. At each iteration k we use the following discretizations of the Dirichlet transmission conditions on each interface

d(l) (d = 1, · · · , 4): P(l),k = P(3−l),k−1 . 3.

(22)

Physical parametrization

The forcing terms appearing in the right-hand sides of Eqs (1) and (2) are specified or parametrized by the unified Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

(a)

(b)

(c)

(d)

1917

(e)

Figure 2. Objective evaluation of 120 h forecasts against radiosondes for the GEM global (blue lines) and Yin-Yang (red lines) global forecast systems for (a) U winds in m s−1 , (b) wind speed in m s−1 , (c) geopotential height (dam), (d) temperature (K) and (e) dew point depression (K). The RMSE (solid lines) and bias (dashed lines) are shown. Verification is performed over North America, for a set of 42 integrations over the winter period 15 December 2008–14 February 2009.

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1918

A. Qaddouri and V. Lee

(a)

(b)

(c)

(d)

(e)

Figure 3. Objective evaluation of 120 h forecasts against radiosondes for the GEM global (blue lines) and Yin-Yang (red lines) global forecast systems for (a) U winds in m s−1 , (b) wind speed in m s−1 , (c) geopotential height (dam), (d) temperature (K) and (e) dew point depression (K). The RMSE (solid lines) and bias (dashed lines) are shown. Verification is performed over North America using KF trigger=0.05, for a set of 42 integrations over the summer period 18 June 2008–18 August 2008.

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

(a)

(b)

(c)

(d)

1919

(e)

Figure 4. Objective evaluation of 120 h forecasts against radiosondes for the GEM global (blue lines) and Yin-Yang (red lines) global forecast systems for (a) U winds in m s−1 , (b) wind speed in m s−1 , (c) geopotential height (dam), (d) temperature (K) and (e) dew point depression (K). The RMSE (solid lines) and bias (dashed lines) are shown. Verification is performed over North America using KF trigger=0.0475, for a set of 42 integrations over the summer period 18 June 2008–18 August 2008.

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1920

A. Qaddouri and V. Lee

(a)

(b)

(c)

(d)

Figure 5. 24 h accumulated precipitation for integration using analysis on 11 August 2010 at 0000 UTC. (a) GEM global (uniform grid), (b) Yin-Yang, (c) Yin-Yang (KF trigger=0.0475) and (d) GEM global (variable grid) with Yin grid in core, all using KF trigger=0.05 except for (c).

(9) the ozone climatology based on ozonesonde and RPN (Recherche Prevision Numerique) physics package. The choice of parametrization depends on the space- and satellite measurements from Fortuin and Kelder time-scales of the forecast. In this project, the physical (Fortuin and Kelder, 1998). processes are the same as in the GEM global forecast system even though we are coupling two LAM grids (the two panels 4. Digital filter and diffusion of the Yin-Yang grid). The parametrization scheme used for the physics configuration is as follows: For both the GEM global and Yin-Yang runs, the diabatic (1) the Interaction Soil Biosphere Atmosphere (ISBA) digital filter (Fillion et al., 1995) with a 6 hour span land-surface scheme for the surface-layer effects was applied at the beginning of the model integration to dampen high-frequency noise from the initial conditions. (B´elair et al., 2003a,b); (2) Geleyn boundary-layer cloud scheme (Geleyn, 1987); In both models, horizontal diffusion is expressed in the (3) Kuo transient shallow convection scheme (Kuo, form of a scale-selective hyper-Laplacian (∇ 6 ) applied to 1974); momentum variables and temperature. This hyper-diffusion (4) the Kain–Fritsch deep convection scheme (Kain and is solved by an implicit scheme described in Qaddouri and Fritsch, 1993); Lee (2008). (5) the Bougeault–Lacarr`ere turbulent mixing-length An enhanced diffusion is applied in the sponge layer scheme (Bougeault and Lacarr`ere, 1989); covering the six uppermost levels. In the GEM global run, (6) the radiative transfer scheme from Li and Barker (Li it uses an implicit horizontal (∇ 2 ) scheme, whereas for the and Barker, 2005); (7) the non-orographic gravity-wave drag scheme by Yin-Yang runs only the explicit nine-point filter scheme was available under the GEM–LAM mode. The coefficient in the Hines (Hines, 1997a,b); (8) the inclusion of a methane oxidation parametrization latter scheme was adjusted to give the closest behaviour to scheme (same scheme used in the ECMWF model); the implicit scheme. c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

1921

(a)

(b)

Figure 6. Quantitative precipitation forecasts for 24–48 h for the GEM global (blue lines with crosses) and Yin-Yang (red lines with stars) global forecasts. The bias is calculated in such a way that when the value is closer to 1.0 there is less bias. The higher threat score indicates more precision in the forecasts. Verification is carried out over North America using the observational surface network from the Standard Hydrometeorological Exchange Format(SHEF) for a set of 42 integrations over the winter period 15 December 2008–14 February 2009.

5. Experimental set-up and evaluation

5.1. Objective evaluation

The goal in this preliminary objective evaluation is to demonstrate that we could use a Yin-Yang grid in the GEM model to produce the same forecast results as the GEM global configuration. The latter uses a global latitude–longitude grid. In order to make a fair comparison in terms of having the closest horizontal resolution, the GEM global had to be configured with twice the number of points along the longitude versus the number of points along the latitude. Our current 33 km operational GEM global has a horizontal grid of 800 × 600 points with 80 vertical levels. We decided to compare a global set-up using a horizontal grid of 800 × 400 points and 80 vertical levels with the Yin-Yang grid system of (600 × 200) × 2 points and 80 vertical levels. Geophysical fields were generated for each respective grid: the GEM global using 800 × 400 grid points and 600 × 200 grid points on each Yin and Yang grid.

For each system (Yin-Yang and GEM global), we conducted two series of tests: 42 winter cases and 42 summer cases where each case is a five-day forecast verified against upper-air and radiosonde observations. Each season spans over a two-month period at every 36 hours. The values of the quantitative precipitation forecasts (QPFs) were also compared between the two systems. This is a typical method used to evaluate models at CMC. When we first compared the objective evaluation of both systems for the winter series against upper-air and radiosonde observations, both gave very similar results, as shown in Figure 2. The blue lines represent the GEM global runs, which are our control runs, and the red lines represent the Yin-Yang runs. Note that the bias and the root-meansquare errors (RMSE) are almost equivalent for both the GEM global runs and the Yin-Yang runs. Each result from each system is almost superimposed for the winds (UU, UV),

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1922

A. Qaddouri and V. Lee

(a)

(b)

Figure 7. Quantitative precipitation forecasts for 24–48 h for the GEM global (blue lines with crosses) and Yin-Yang (red lines with stars) global forecast. The bias is calculated in such a way that when the value is closer to 1.0 there is less bias. The higher threat score indicates more precision in the forecasts. Verification is carried out over North America using the observational surface network from the Standard Hydrometeorological Exchange Format (SHEF) for a set of 42 integrations over the summer period 18 June 2008–18 August 2008.

temperature (TT), geopotential height (GZ) and dewpoint depression (ES). We then ran the summer series, where the results in Figure 3 showed a slight unfavourable significant difference for the Yin-Yang system from the bias in the GZ above 100 mb. On a closer inspection, we noticed that the YinYang GEM model was underpredicting larger amounts of precipitation, particularly for rainfall over 20 mm in a 24 h period and specifically over land. It was also noted that the vertical motion was a bit less intense (smoother) in the Yin-Yang than in the GEM global configuration. It is known that the intensity of the vertical motion is directly related to the resolution of the grid. The factor in the Kain–Fritsch (KF) deep convection scheme is usually adjusted to suit the average resolution of the grid in order to control how easily the KF scheme can be triggered. The finer the resolution, the higher the grid-scale vertical velocity. When the temperature of a buoyant parcel plus temperature perturbation associated with the grid-scale vertical velocity c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

is higher than the environmental temperature, then the KF scheme is triggered. Unfortunately, in our current operational models the KF trigger value is tuned to be dependent on an averaged horizontal resolution and not on the local resolution nor on the true latitude. We found that lowering the KF trigger value slightly to 0.0475 m s−1 was sufficient to obtain the same validation as the GEM global for the summer series (Figure 4). Note that there is no longer a significant difference in the bias of the GZ, as well as there being a slight improvement in the bias of TT (temperature) between 400 and 250 mb. The winter series for the Yin-Yang model were repeated with the newly adjusted KF trigger value of 0.0475 m s−1 and this had no impact on the results. One of the biggest challenges in comparing the two systems is to match the grids as closely as possible. The Yin-Yang grid remains quasi-uniform while the global grid has an increasing resolution for grid points approaching the poles. The grids in GEM are designed in such a way that Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

1923

(a)

(b)

Figure 8. The 72 h forecast of the air temperature (TT) close to the surface for (a) GEM global (uniform grid) and (b) Yin-Yang grid.

even with almost the same resolution at the Equator, the GEM global grid points are staggered to the grid points of the Yin-Yang grid. To test the sensitivity of precipitation in the location of grid points, we decided to simulate a global run using the Yin grid as the core domain in a global variable grid described in Cˆot´e et al. (1993). Figure 5 shows four plots of the same region of the USA with a 24 h forecast snapshot of accumulated precipitation. All use a KF trigger of 0.05 m s−1 except for Figure 5(c). Note that the precipitation maxima of the global uniform grid shown in Figure 5(a) are much greater than those in the original Yin-Yang run shown in Figure 5(b). In Figure 5(c), the maxima of the precipitation increased for the Yin-Yang grid, due to decreasing the KF trigger to 0.0475 m s−1 . In Figure 5(d), the precipitation shown is from a run using the global variable grid using the Yin grid as its core uniform domain and it also used the original KF trigger of 0.05 m s−1 . Note that Figure 5(d) bore the greatest resemblance to the original Yin-Yang run in Figure 5(b) because the grid points were the most closely matched. This indicates that the precipitation pattern from the GEM can be affected not only by the horizontal resolution of grid points but also by their location. One of the ways of validating precipitation at CMC is to look at the quantitative precipitation forecasts (QPFs) from each model of both the winter cases (Figure 6) and the summer cases (Figure 7). In each figure there are two c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

plots: one is the bias and the other is the threat score, where the blue represents GEM global and the red represents YinYang. The bias calculated is based on the number of forecasts for a category divided by the number of observations in this category. The threat score is a measure of the number of correct forecasts for a category divided by the number of events forecast or observed in this category. In looking at these two QPF plots, it seems that there are no significant meteorological differences. A major advantage in using the Yin-Yang grid system to represent the Earth’s spherical geometry is that one avoids singularities and mesh convergence in the polar regions. There is an expectation that a Yin-Yang model may give better results in the regions near the poles. In all the integrations done with the GEM global model, we can find numerical noise in the pole areas and an example is shown in Figure 8(a) for the 72 h forecast of air temperature (TT) close to the surface. The same forecast produced by the YinYang model shows much smoother contours at the south pole and maintains the same pattern as the GEM global elsewhere (see Figure 8(b)). 5.2. Preliminary parallel performance As in the experimental set-up, we have set the resolution to be the same at the Equator for both systems (Yin-Yang and GEM global). We have chosen to use fast Fourier transforms Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1924

A. Qaddouri and V. Lee

Figure 9. Wall-clock time in seconds and speed-up versus number of CPUs for 96 time steps. Both GEM global and Yin-Yang runs are configured for 39 km resolution at the Equator.

(FFT) in the EBV solver for both systems, which meant that the choice of the number of longitude points (N) is limited. In the case of Yin-Yang, there is an additional constraint whereby the number of longitude points (N) must be three times the number of latitude points (M). For the global grid, N is twice M. The resolution chosen for this preliminary parallel performance test is 39 km where the global grid is 1024 × 512 and each Yin-Yang subgrid is 799 × 267 points. Both have the same definition for the 80 vertical levels and both use the same configuration for the physics. The time step given is 720 s. For each test, we let the model run 96 time steps and measured the wall-clock time. We tried various processor topologies for the two systems and retained the best ones for each system. Our best approach to compare performance, especially in the context of operational runs, is to plot the wall-clock time versus the number of CPUs used in each run. In Figure 9(a), we see plots of both systems with the wall-clock time of a run versus the number of CPUs used in a given run. The red plot line represents the Yin-Yang runs, and this starts to outperform the GEM global run (blue line) when using more than 512 CPUs. Another way to compare the parallel scalability of the two models is to plot the parallel speed-up versus the number of CPUs used. In Figure 9(b), the speed-up of 1.0 is based on the runs c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

using 192 CPUs, as we cannot run a 39 km resolution model on 1 CPU. The results plotted are the best timings for the number of CPUs given to each system, since it is rare that both can use the same number of CPUs. This figure shows a better overall scalability with the Yin-Yang model. Looking at Table 1, we can see that the advection scales better in the Yin-Yang case because of the LAM mode when we are increasing the number of processors. At the same time, the smaller value in the physics computation indicates that there are fewer points to compute for the physical processes. The disadvantage of Yin-Yang is in the computation for the solver, as it must solve iteratively between Yin and Yang to obtain its solution. The other additional cost in the YinYang system is the extra communication between the two panels in the update of the boundary conditions and other intermediate steps to complete the two-way coupling. The code still needs further optimization and we hope to see an even better speed-up in the future. 6.

Conclusion

The future global forecasting Yin-Yang model has been developed by a two-way coupling method between two GEM LAMs discretized on the two panels of the Yin-Yang grid. A Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

Canadian GEM model on the Yin-Yang grid

1925

Table I. Performance breakdown of timings in seconds. Mode Yin-Yang Global Yin-Yang Global Yin-Yang Global Yin-Yang Global

PE Topo

NEST

BAC

SOL

ADV

PHY

TOT

CPUs

(12 × 8 × 1) × 2 (3 × 32 × 2) (20 × 8 × 1) × 2 (2 × 40 × 4) (16 × 12 × 1) × 2 (3 × 32 × 4) (16 × 8 × 2) × 2 (4 × 37 × 4)

64 n/a 54 n/a 43 n/a 32 n/a

38 6 19 5 19 4 16 4

122 39 125 35 73 28 67 23

134 162 93 127 74 109 65 82

233 272 130 222 113 145 90 115

676 597 496 508 386 388 326 325

192 192 320 320 384 384 512 592

PE Topo –processor topology (Npex × Npey × OpenMP). NEST –exchange of boundary conditions. BAC –Back Substitution. SOL –Solver. ADV –Advection. PHY –Physics.TOT –Total Wallclock. CPUs –Total Number of CPUs used.

comparison has been made by matching as closely as possible the grid resolution of each system (Yin-Yang and GEM global) and then performing a model validation for both systems. Preliminary results of the new Yin-Yang system show that it is capable of producing meteorological forecasts equivalent to the GEM global model. One big advantage in using this type of system is to avoid the singularities and convergence of meridians in the polar regions. Another benefit is to have a more balanced computational load for scalability purposes. Our next step would be to compare objective evaluations with a much higher resolution grid and using a non-hydrostatic version of the GEM model. The results for parallel performance also indicate better scalability compared with the global model when the system is given more processors. Investigation and development will continue, to consider whether this Yin-Yang system could become the future global forecasting system at CMC. Acknowledgement ˆ e and Claude The authors thank their colleagues Jean Cot´ Girard and the two anonymous reviewers for their helpful comments. References Adcroft A, Campin J-M, Hill C, Marshal J. 2004. Implementation of atmosphere–ocean general circulation model on the spherical cube. Mon. Weather Rev. 132: 2845–2863. Baba Y, Takahashi K, Sugimura T. 2010. Dynamical core of an atmospheric general circulation model, on Yin-Yang grid. Mon. Weather Rev. 138(10): 3988–4005, DOI:10.1175/2010MWR3375.1. B´elair S, Crevier L-P, Mailhot J, Bilodeau B, Delage Y. 2003a. Operational implementation of the ISBA land surface scheme in the Canadian regional weather forecast model. Part I: Warm season results. J. Hydrometeorol. 4: 352–370. B´elair S, Brown R, Mailhot J, Bilodeau B, Crevier L-P. 2003b. Operational implementation of the ISBA land surface scheme in the Canadian regional weather forecast model. Part II: Cold season results. J. Hydrometeorol. 4: 371–386. B´elair S, Roch M, Leduc A-M, Vaillancourt P-A, Laroche S, Mailhot J. 2009. Medium-range quantitative precipitation forecasts from Canada’s new 33 km deterministic global operational system. Weather and Forecasting 24: 690–708. Bougeault P, Lacarr`ere P. 1989. Parameterization of orography-induced turbulence in a mesobeta-scale model. Mon. Weather Rev. 117: 1872–1890. Cˆot´e J. 1988. A Lagrange multiplier approach for the metric terms of semi-Lagrangian models on the sphere. Q. J. R. Meteorol. Soc. 114: 1347–1352. Cˆot´e J, Staniforth A. 1988. A two-time-level semi-Lagrangian semiimplicit scheme for spectral models. Mon. Weather Rev. 116: 2003–2012. c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Cˆot´e J, Roch M, Staniforth A, Fillion L. 1993. A variable-resolution semi-Lagrangian finite-element global model of the shallow-water equations. Mon. Weather Rev. 121: 231–243. Cˆot´e J, Gravel S, M´ethot A, Patoine A, Roch M, Staniforth A. 1998. The operational CMC–MRB Global Environmental Multiscale (GEM) model. Part I: Design considerations and formulation. Mon. Weather Rev. 126: 1373–1395. Dudhia J, Bresch J-F. 2002. A global version of PSU-NCAR mesoscale model. Mon. Weather Rev. 130: 2989–3007. Fillion L, Mitchel H-L, Ritchie H-R, Staniforth A. 1995. The impact of a digital filter finalization technique in a global data assimilation system. Tellus 47A: 304–323. Fortuin P, Kelder H. 1998. An ozone climatology based on ozonesonde and satellite measurements. J. Geophys. Res. 103: 31709–31734. Geleyn J-F. 1987. ‘Use of a modified Richardson number for parameterizing the effect of shallow convection. Short- and medium-range numerical weather prediction’. In WMO/IUGG NWP Symposium. WMO: Tokyo, Japan; pp 141–149. Girard C, Plante A, Gravel S, Qaddouri A, Chamberland A, Spacek L, Lee A, Desgagn´e M. 2010. GEM4.1: A non-hydrostatic atmospheric model (Euler equations). RPN document. http://collaboration. cmc.ec.gc.ca/science/rpn/publications/pdf/GEM4.1.pdf. Hines C-O. 1997a. Doppler-spread parameterization of gravity-wave momentum deposition in the middle atmosphere. Part 1: Basic formulation. J. Atmos. Solar–Terr. Phys. 59: 371–386. Hines C-O. 1997b. Doppler-spread parameterization of gravity-wave momentum deposition in the middle atmosphere. Part 2: Broad and quasi monochromatic spectra, and implementation. J. Atmos. Solar–Terr. Phys. 59: 387–400. Kageyama A, Sato T. 2004. The ‘Yin-Yang grid’: An overset grid in spherical geometry. Geochem. Geophys. Geosystems 5: Q09005. DOI:10.1029/2004GC000734. Kain J, Fritsch J-M. 1993. Convective parameterization for mesoscale models: The Kain–Fritsch scheme. The representation of cumulus convection in numerical models. Meteorol. Monogr. 46: 165–170. Kuo H-L. 1974. Further studies of the parameterization of the influence of cumulus convection on large-scale flow. J. Atmos. Sci. 31: 1232–1240. Li J, Barker H-W. 2005. A radiation algorithm with correlated-k distribution. Part I: Local thermal equilibrium. J. Atmos. Sci. 62: 286–309. Qaddouri A. 2008. Optimized Schwarz methods with the Ying-Yang grid for shallow water equations. In Domain Decomposition Methods in Science and Engineering XVII. Lecture Notes in Computational Science and Engineering 60 347–353. Qaddouri A. 2011. Nonlinear shallow-water equations on the Yin-Yang grid. Q. J. R. Meteorol. Soc. 137: 810–818, DOI:10.1002/qj.792. Qaddouri A, Lee V. 2008. ‘Solution of the implicit formulation of high order diffusion for the Canadian Atmospheric GEM model’. In Proc. 2008 Spring Simulation Multiconf., High Performance Computing Symp., Hamilton JA Jr et al. (eds). Soc. For Modeling and Simulation Internat.: Ottawa, Canada; pp 362–367. Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)

1926

A. Qaddouri and V. Lee

Qaddouri A, Laayouni L, Loisel S, Cˆot´e J, Gander M-J. 2008. Optimized Schwarz methods with an overset grid for the shallowwater equations: Preliminary results. Appl. Numerical Math. 58: 459–471. Rivest C, Staniforth A, Robert A. 1994. Spurious resonant response of semi-Lagrangian discretizations to orographic forcing: Diagnosis and solution. Mon. Weather Rev. 122: 366–376.

c 2011 Crown in the right of Canada. Published by John Wiley & Sons Ltd.

Yanenko N-N. 1971. The method of Fractional Steps. Springer: Berlin; 160 pp. Yeh K-S, Cˆot´e J, Gravel S, M´ethot A, Patoine A, Roch M, Staniforth A. 2002. The CMC–MRB global environmental multiscale (GEM) model: Part III – Nonhydrostatic formulation. Mon. Weather Rev. 130: 339–356.

Q. J. R. Meteorol. Soc. 137: 1913–1926 (2011)