New Version of the CMA-GFS Dynamical Core Based on the Predictor–Corrector Time Integration Scheme

Numerical weather prediction (NWP) has become the key tool in modern weather forecasting. NWP is the technique used to obtain an objective forecast of the weather up to possibly two weeks into the future by solving the governing equations that describe the evolution of variables that define the present state of the atmosphere. Since 2001, the China Meteorological Administration (CMA) has been independently developing a new-generation NWP system, which is called the Global/Regional Assimilation and Prediction System (GRAPES) (Xue and Liu, 2007; Xue and Chen, 2008).

The grid-point-based GRAPES model uses a regular latitude–longitude grid in the horizontal and a terrain-following height-based grid in the vertical, C-grid staggering in the horizontal and Charney–Phillips staggering in the vertical, and employs a classical two-time-level (2TL) semi-implicit semi-Lagrangian (SISL) scheme for time discretization (Xue and Chen, 2008). To date, a complete operational numerical prediction system has been developed by the CMA, based on the GRAPES model, which includes the CMA global forecast system (CMA-GFS) global 10-day deterministic forecast (25 km) (Shen et al., 2017), global 15-day ensemble forecast (CMA-GEPS, 50 km) (Chen and Li, 2020), China 36-h high-resolution forecast (CMA-MESO, 3 km) (Huang et al., 2017), China 3-day regional ensemble prediction system (CMA-REPS, 10 km) (Wang et al., 2018), and Indian Ocean–Northwest Pacific 5-day typhoon forecast (CMA-TYM, 9 km) (Zhang et al., 2017).

The semi-Lagrangian (SL) technique is unconditionally stable, but its combination with the semi-implicit (SI) scheme succeeds in maintaining the stability of a fully coupled system of atmospheric equations. Consequently, the SISL approach has been used widely in operational NWP for more than 30 years (Beljadid et al., 2014). The 2TL SISL schemes are notably advantageous for large-scale NWP because they remain stable when long time steps are used, and are twice as efficient in comparison with three-time-level schemes (Gravel et al., 1993). However, as identified by Gravel et al. (1993), several issues must be handled with special attention when implementing a 2TL SISL scheme, e.g., the selection of the reference state in the linearization, and the discretized formulation of advection and nonlinear terms (Ritchie and Tanguay, 1996; Simmons and Temperton, 1997; Temperton et al., 2001). It should be noted that these issues are not obvious in a three-time-level SISL scheme, possibly owing to the application of the Asselin filter.

Usually, implementation of the SI algorithm requires a certain reference state to divide the nonlinear terms into linear terms and nonlinear residuals. The implicit treatment applies to the linear terms, while the nonlinear residuals are traditionally handled by temporal extrapolation. However, this extrapolation has been proven to induce instability in some circumstances (Gospodinov et al., 2001; Hortal, 2002; Bénard, 2003; Durran and Reinecke, 2004).

Several operational NWP centers have been successful in using SL advection (Bates et al., 1993). As discussed by both Davies et al. (2005) and Cordero et al. (2005), the traditional treatment of the SL trajectory calculation adopts the midpoint way, which requires temporal extrapolation of the midpoint winds at half a time interval in the future, which might also cause instability (McDonald, 1997; Hortal, 1998). Furthermore, as shown by Cordero et al. (2005), temporal extrapolation of the midpoint winds at $ {t}^{n+\frac{1}{2}} $ might substantially distort the vertical structure of the normal modes in a nonhydrostatic SISL model.

In a brief summary, sources of instability in the classical 2TL SISL scheme have been identified from two temporal extrapolations related to the SI manipulation of nonlinear residuals and the calculation of midpoint winds of the SL trajectory. As commented by Bénard (2003), the classical 2TL SISL scheme might not be suitable for the fully compressible equations, and if high-level accuracy is desired with a 2TL SISL scheme at high resolution, a more robust approach must be used. This prompted the authors to seek more a robust scheme for improving the CMA-GFS forecast skill to meet the needs of high-resolution operational applications.

Hortal (2002) proposed the Stable Extrapolation Two-Time-Level Scheme as a contribution to the solution of the instability problem. However, Durran and Reinecke (2004) showed that such a scheme was not a perfect choice because, except at the origin, the imaginary axis lies outside the region of absolute stability. Bénard (2003) showed that iteration of the SI discrete set of equations could be an appropriate approach to enhance stability, and this iterative approach has been employed in the Canadian operational Global Environmental Multiscale model and the Unified Model of the U.K. Met Office (Diamantakis et al., 2007).

Several studies have proposed a number of explicit predictor–corrector schemes for use in atmospheric models (Kurihara, 1965; Durran and Blossey, 2012; Kar, 2012), in which the nonlinear residual terms are handled by a multistage approach. In this paper, we seek improved time discretization by considering predictor–corrector schemes with an implicit method at each stage, and by combining the SI predictor–corrector scheme with SL advection. Additionally, appropriate choice of the reference profile in constructing the SI algorithm is important for accuracy and convergence of the SISL scheme (Simmons et al., 1978; Simmons and Temperton, 1997; Temperton et al., 2001; Bénard, 2003). Usually, the reference profile is considered stationary in time. Fundamentally, the nonlinear residuals should be small when choosing the reference profile. In the hydrostatic primitive equation context, Simmons and Temperton (1997) showed that 2TL SISL schemes have stability constraints that are more stringent than those of their three-time-level counterparts in choosing the reference profile. In the nonhydrostatic equation context, the situation might become even more complicated. In the nonhydrostatic Unified Model of the U.K. Met Office, a dynamically determined reference profile is used (Diamantakis et al., 2007). Therefore, selection of an appropriate SI reference profile that can improve the CMA-GFS dynamics is revisited in this paper.

The remainder of this paper is structured as follows. In Section 2, an overview of the dynamic research progress of the CMA-GFS is given, which includes description of both the classical 2TL SISL scheme and the new predictor–corrector algorithm, and discussion on the new implementation of the 3D reference profile and the hybrid vertical coordinate system. In Section 3, the modified discrete model equations are presented. Results from idealized tests and retrospective forecasts using the CMA-GFS are presented in Section 4. Finally, the derived conclusions are summarized in Section 5.

First, a brief introduction is given to the governing equations and the current temporal integration scheme of the CMA-GFS. Second, details of the components of the new research progress are introduced. Finally, a summary of the version upgrade content is provided.

The prognostic variables of the CMA-GFS are the 3D wind components $ \left(u,v,w\right) $, potential temperature $ \left(\theta \right) $, Exner pressure $ \left(\Pi \right) $, and mixing ratios of the moisture equations (q_X, X: various phases of water), as set out in Xue and Chen (2008). The equations, which are nonhydrostatic with the shallow atmospheric approximation in spherical polar coordinates $ \left(\lambda ,\phi ,r\right) $ with the origin at the earth’s center, can be expressed as follows:

where $ {\theta }_{v} $ is virtual potential temperature, $ f $ is the Coriolis force constant, $ g $ is gravitational acceleration, $ {C}_{P} $ is specific heat capacity at constant pressure, $ \gamma ={C}_{P}{R}^{-1} $, $ R $ is the gas constant, $ {\delta }_{NH} $ is the switch for the hydrostatic balance, $ {D}_{3} $ is 3D divergence, and $ {F}_{\theta }^{*} $ is the source and sink of heat.

In implementing the SI scheme, the prognostic equations in Section 2.1 are linearized around a stationary reference state. The linearized prognostic equations can be written in the following general form:

where $ X $ denotes position, $ \psi $ represents any of the prognostic variables, $ L $ represents the linear terms, and $ N $ represents the nonlinear terms. The first version of the GRAPES dynamical core utilized the classical 2TL SISL scheme, which can be written as follows:

where $ n $ is the time level, $ \Delta t $ is the time step, $ d $ denotes the departure point, and $ \alpha $ is a time-weighted coefficient.

The classical 2TL SISL scheme needs to handle the temporal extrapolation twice to estimate the nonlinear residuals $ N $ at the $ n+1 $ time step and the midpoint winds at the $ n+\dfrac{1}{2} $ time step to obtain the SL trajectory (McDonald and Bates, 1987). In the first version, the following formulas are adopted for $ {N}^{n+1} $ and $ V(X,{t}^{n+\frac{1}{2}}) $, while $ V(X,{t}^{n+\frac{1}{2}}) $ is used to temporally interpolate to the mid-point $ \dfrac{{X}_{a}+{X}_{d}}{2} $ of the trajectory to obtain the wind fields:

As mentioned above, the 2TL extrapolations have been identified as sources of instability. Meanwhile, the SL technique suffers from spurious SL orographic resonance (Tanguay et al., 1992; Temperton et al., 2001). Tanguay et al. (1992), Rivest et al. (1994), and Ritchie and Tanguay (1996) proposed and examined a technique of spatial averaging along the SL trajectory to reduce such noise. Therefore, to maintain stability in the classical 2TL SISL scheme, $ \alpha $ usually takes a value larger than $ \dfrac{1}{2} $, which damps the scheme to first-order accuracy temporally. The use of a value of $ \alpha $ larger than $ \dfrac{1}{2} $ induces dramatic loss in kinetic energy with increasing forecast range, and damps the transient features of weather systems.

To improve the stability of the CMA-GFS, an SI predictor–corrector scheme is proposed. As described in Section 2.2, a general Lagrangian-form differential equation is considered here, and the source and sink terms are ignored. For clarity, we separately discuss the schemes for handling the linear and nonlinear residual terms first, and then extend its use to 2TL SISL discretization.

For the nonlinear residual terms, many possible predictor–corrector discretization schemes that avoid the temporal extrapolation are available, e.g., the Adams–Bashforth trapezoidal method and the Adams–Bashforth–Moulton method (Durran and Blossey, 2012; Kar, 2012). These schemes, which have been discussed widely in textbooks or in relation to application to the numerical solution of ordinary differential equations, are usually called linear multistep methods.

To retain the currently implemented SISL scheme in the new algorithm, the Adams–Bashforth predictor followed by a trapezoidal corrector is adopted as follows:

which has second-order accuracy.

For the linear terms, there are many implicit discretization methods available, as discussed in Durran and Blossey (2012). Adopting a stable and accurate scheme for the linear terms is crucial, especially in relation to high-resolution NWP. Moreover, such a scheme should be compatible with the scheme proposed above for the nonlinear residual terms. To fully exploit the original SI implementation in the CMA-GFS, we use a modified second-order Adams–Moulton formulation suggested by Durran and Blossey (2012) in the first step and trapezoidal averaging in the second step:

The nature of the stability of the above second-order approximation of the linear terms has been discussed by Durran and Blossey (2012), Clancy and Pudykiewicz (2013), and Beljadid et al. (2014).

The combination of the above explicit predictor–corrector scheme with a two-step implicit discretization constitutes the SI predictor–corrector scheme of this work. To apply the above SI predictor–corrector scheme into the 2TL SISL framework, the prognostic equations may be written symbolically as follows:

where $N{\left({\psi }^{*}\right)}^{n+\frac{1}{2}}=\dfrac{3}{2}N\left({\psi }^{n}\right)-\dfrac{1}{2}N\left({\psi }^{n-1}\right)$, and generally ${\alpha }_{1}=\dfrac{3}{4}\;{\rm{and}} \hspace{0.33em}{\alpha }_{2}=\dfrac{1}{2}$.

It is obvious that the above formulae are very similar to the original handling of the 2TL SISL scheme in the CMA-GFS, and that we can easily extend the original scheme into the more stable and accurate SI predictor–corrector scheme.

In the SI predictor–corrector version of the CMA-GFS, which is presented in this paper, the departure point is recomputed at the corrector stage. Using the midpoint approximation, the departure point is estimated at the predictor stage as follows:

where the velocity field at time $ {t}^{n}+\dfrac{\Delta t}{2} $ is estimated by the temporal extrapolation, as described in Section 2.2. The tag at the top of the vector is omitted.

At the corrector stage, the departure point is evaluated by using the following formula:

where ${V}^{*}\left(\dfrac{{x}_{a}+{x}_{d}^{*}}{2}\right)$ is the spatial interpolation to the midpoint of the wind vector at the predictor stage, and ${V}^{n}\left(\dfrac{{x}_{a}+{x}_{d}}{2}\right)$ is the spatial interpolation to the midpoint of the wind vector at the previous time step.

As mentioned above, an appropriate choice of the reference profile for implementing the SI algorithm is important for accuracy and convergence of the SISL scheme. Simmons et al. (1978) and Simmons and Temperton (1997) examined the stability in terms of the atmospheric and reference isothermal temperature profiles in the hydrostatic equation context. Subsequently, Bénard (2003) further generalized the issues for any usual meteorological system of NWP interest including the nonhydrostatic equations. It should be noted that computational instability might arise with an SI scheme if there is marked difference between the actual temperature profile and the reference profile. Ideally, the nonlinear residuals of the SI scheme should be as small as possible, such that the linear part of the SI scheme can remain close to the tangent–linear operator of the governing equations.

The first version of the CMA-GFS adopted the isothermal stationary atmosphere as the reference state. However, this might induce reasonably large nonlinear residuals, especially in the upper levels, and it might influence computational accuracy and cause instability.

Su et al. (2018, 2020) investigated the implementation of a 3D reference profile into the CMA-GFS SI discretization. The 3D reference profile can be derived from the climatology or the initial values using the hydrostatic balance relationship between $\tilde \theta $ and $\tilde \Pi$. In some cases, $ \tilde{\theta } $ must be adjusted to ensure that it increases monotonically in the vertical direction. In comparison with the 1D reference, the right-hand side of the discretized equations will include extra terms that take account of the horizontal variations of the reference state.

Utilization of a terrain-following coordinate system is the standard choice in NWP and climate modeling for simple treatment of the lower boundary conditions (Husain et al., 2019). However, gradual flattening of the coordinate surfaces in the upper levels of a model is important for reducing errors in the calculations of the pressure gradient force and advection transport (Simmons and Burridge, 1981). Usually, a hybrid approach is adopted in terrain-following coordinate systems. In the first version of the CMA-GFS, a height-based terrain-following coordinate system was utilized but without flattening of the coordinate surfaces in the upper levels (Gal-Chen and Somerville, 1975). Recently, a height-based terrain-following hybrid vertical coordinate system was introduced to accompany implementation of the SI predictor–corrector scheme:

where $ z $ is geopotential height, $ \widehat{z} $ represents the coordinates of the calculation domain after vertical transformation, $ {Z}_{s} $ is the height of orography, $ {Z}_{T} $ is the height of the model top, a is the coefficient for specifying and adjusting the flattening levels, $ {z}_{c1} $ is the height at which the hybrid coordinate begins, and $ {z}_{c2} $ is the height at which the coordinate begins to flatten.

Equation (13) represents the generic representation of the predictor–corrector SISL discretization applied in the new CMA-GFS model. For clarity, we further describe the equations at the predictor and corrector stages to clearly show the differences from the original SISL scheme. This section also demonstrates that only a few simple changes to the model code are necessary for implementing the predictor–corrector scheme. Again, we describe only the adiabatic part of the scheme because the physics parameterizations and physics/dynamics coupling algorithm do not make any direct changes to the updated dynamical core.

In the predictor step, the following generic equation can be applied to each discrete equation for the prognostic variables, and it retains the same form as the original presented in Xue and Chen (2008). The resulting scheme may be written as follows:

In comparison with the original formulas, i.e., (2.3.12), (2.3.14), (2.3.51)–(2.3.53) in Xue and Chen (2008), the only difference is in the second term of the right-hand side of Eq. (17), which can be simply modified in the code. Additionally, it should be noted that the linear term L should include extra terms from the 3D reference profile. In this step, a linear Helmholtz equation must be solved, but the solver costs can be controlled only if the solution could be used as a good estimate for the corrector step. Moreover, a much cheaper scheme for computation of the advection of the moisture variables could be applied in this step, e.g., the quasi-monotonic SL scheme.

In the corrector step, the resulting scheme may be written as follows:

Again, the only difference from the original scheme is in the second term of the right-hand side of Eq. (18). It should be noted that the temporal weighting (or off-centered parameter) $ {\alpha }_{2} $ takes a value of 0.5 in the corrector step to achieve second-order accuracy. Finally, a linear Helmholtz equation is solved, and the solver costs should be smaller because the initial guess for the Helmholtz solver should be closer to the solution, meaning that the solver should converge faster.

In addition to the major improvements mentioned above, some minor improvements have also been implemented. 1) With the new spectral filter terrain, better distribution of the kinetic energy spectrum is obtained on the shortwave scale, which provides necessary information for further improvements in small-scale circulations and terrain-forced precipitation. 2) With optimization of artificial restrictions, excessive damping and smoothing are removed, which improves simulation capabilities for small-scale systems.

Integration of the above improvements forms the new version of the dynamical core, which combined with the corresponding modifications of the 4DVAR system, represents version 3.0 of the CMA-GFS. The CMA-GFS 3.0, which has been operational since June 2020, replaced the previous version of the model (i.e., CMA-GFS 2.4). Details of the differences between the two versions are listed in Table 1. In this version upgrade, the physics process was not changed.

This section summarizes the added value of the new version of the dynamical core in comparison with the original in terms of accuracy, stability, and efficiency.

The main purpose of the Dynamical Core Model Intercomparison Project (DCMIP) is to compare the performance of different dynamical cores through a series of benchmark tests (Ullrich et al., 2013), among which test No. 2_0, assessing the accuracy of the PGF over steep terrain in a stationary atmosphere, is considered in this work. A moderately steep Schar-like circular mountain was set in the tropics with its center at 0°, 150°W. The radius and maximum height of the mountain were 135° and 2000 m, respectively (Fig. 1). The horizontal resolution was 1.0° and 30 vertical layers with uniform depth of 400 m were used. The simulation was initialized as a steady-state atmosphere and integrated over 6 days with a time step interval of 600 s. Details of the test design and relevant parameters can be found in Ullrich et al. (2013).

From an analytic perspective, the model atmosphere should retain the steady state adopted as the initial condition. However, the coordinate transformation and numerical discretization over steep terrain introduce errors in the calculation of the PGF, and then trigger certain movement. Therefore, this test provides insight into the accuracy of the pressure gradient in models with terrain-following coordinates. The magnitude of the movement could be considered to represent the errors from the calculation of the PGF over the terrain. Comparing the results obtained from the DCMIP test reveals that the error in the PGF of the previous operational version is large, and that the horizontal wind speed (Fig. 2a) is up to 2 m s⁻¹, although both the zonal and the vertical velocity fields (Fig. 2c) show good correlation with the terrain; conversely, the new dynamic version dramatically reduces the error, i.e., the horizontal wind speed (Fig. 2b) is < 0.02 m s⁻¹, comparable with that of the Icosahedral Nonhydrostatic model, which is at the top of the list in terms of performance in this test. The structure of both the zonal and the vertical velocity (Fig. 2d) shows certain correlation with the terrain, but the distribution of each is discontinuous, similar to the test results obtained using the Model for Prediction Across Scales.

Figure 1.  Vertical section along the equator of the terrain used in DCMIP test No. 2_0 (m).

Figure 2.  Vertical section of (a, b) the zonal wind (m s⁻²) and (c, d) the vertical speed (m s⁻²) along the equator in DCMIP test No. 2_0 after integration for 6 days using the (a, c) previous operational version and (b, d) new dynamic version. Note that the color scale in (b, d) is at least one degree less than that in (a, c).

Figure 3b shows the 700-hPa geopotential height after integration over 15 days from the previous operational and new dynamic versions. It can be seen that both versions can replicate the Rossby wave pattern (Fig. 3b) presented in Jablonowski et al. (2008). It is also found that the wave excitation, propagation, and extinction processes were well simulated (not shown). However, in comparison with the previous operational version, the new dynamic version produces a stronger wave downstream of the mountain (Fig. 3b), especially with the 600-s time step, which is a result that is much closer to the reference in Jablonowski et al. (2008).

Figure 3.  (a) Terrain height (m) along 30°N and (b) geopotential height (gpm) at the 700-hPa level after integration over 15 days. Black solid and dashed lines in (b) represent the results from the previous operational version using a time step of 600 and 1200 s, respectively. Red lines are the same as the black lines but for the results from the new dynamic version.

The traditional SL method based on point mapping cannot guarantee mass conservation in a long-term integration process (Su et al., 2016). This is clearly observed in Fig. 4, in which the total atmospheric mass is reduced by approximately 0.1% after 15 days relative to its initial value in the previous operational version. However, the new dynamic version substantially improves model performance in terms of mass conservation, especially when a short time step is used. In real-case NWP, the mass loss during long-term integration in the CMA-GFS is even more severe; therefore, the advantage of the new dynamic version will be even more evident.

Figure 4.  Percentage change in total atmospheric mass relative to initial value during integration. Black solid and dashed lines represent the results from the previous operational version using a time step of 600 and 1200 s, respectively. Red lines are the same as the black lines but for the results from the new dynamic version.

Comparing the KES of a numerical model with that observed is an effective approach for evaluating model dissipation (Zheng et al., 2008). It is well known that the KES in the free troposphere and lower stratosphere is characterized by a power law with slope of k⁻³ for large-scale systems and of k^−5/3 for small-scale systems. To calculate the KES, 12 simulations were conducted, which were initialized at 0000 UTC on the first day of each month in 2008 using the NCEP final analyses (FNL) data. These simulations used a horizontal resolution of 0.25° and a time step of 300 (450) s for the previous operational version (new dynamic version). The mean KES at three isobaric heights (i.e., 100, 500, and 850 hPa) was calculated with the model outputs after 24-h integration from both the previous operational version and the new dynamic version.

As shown in Fig. 5, the KES at the three isobaric heights from the new dynamic version drops less rapidly than that from the previous operational version, especially at 850 hPa, which indicates that energy dissipation is weaker in the new dynamic version. This might be helpful in alleviating the problem of small- and meso-scale systems that were too weak in the previous operational CMA-GFS.

Figure 5.  KES at 100, 500, and 850 hPa. Solid line: previous operational version; long dashed line: new dynamic version; short dashed line: reference for actual atmosphere; dotted line on the right: position of truncated wavenumber 360, which corresponds to 4${\Delta }\mathit{x}$ at 0.25°.

One of the clearest benefits of the new dynamic version is increase in atmospheric variability, which is realized through improvement in accuracy and reduction in artificial restrictions. This can be seen from the globally integrated mass-weighted eddy kinetic energy (MEKE) (Regan et al., 2020), $\mathrm{M}\mathrm{E}\mathrm{K}\mathrm{E}=\sum \dfrac{1}{2}m\left({u'}^{2}+{v'}^{2}\right)$, where $ {u}' $ and ${v}' $ represent the eddy velocity field after deducting the monthly mean field, and m is the mass at each grid point calculated by multiplying the volume by the density. To calculate MEKE, a set of 30 3-day real-case simulations were conducted using the NCEP FNL data. In tests with each set of resolutions, the new dynamic version had a time step 1.5 times that of the previous operational version.

It can be seen from Fig. 6 that the new dynamic version can simulate more abundant eddy motion, and that the difference between the two versions gradually narrows as resolution increases. Similar to the results regarding the KES, this finding indicates that the new dynamic version has stronger capability in simulating eddies and small-scale motion.

Figure 6.  MEKE for various resolutions. Solid line and crosses represent the previous version; long dashed line and circles represent the new dynamic version.

A baroclinic wave is triggered when using steady-state initial conditions with an overlain zonal wind perturbation, and the numerical model should simulate the occurrence, development, explosion, and break of baroclinic waves (Jablonowski et al., 2008). For this test, the model resolution was 0.5° and 60 stretched vertical layers were used. Three time steps were tested, i.e., 600, 1200, and 1800 s, and the integration period was 15 days.

Figure 7 shows the temperature at the 850-hPa level after integration over 9 days from the previous operational and new dynamic versions. It can be seen that both versions simulate the distribution of baroclinic waves reasonably well in comparison with the results in Jablonowski et al. (2008), but there are also obvious differences: the intensity of the fluctuations simulated by the previous operational version cannot be maintained when the time step is doubled and tripled, and obvious attenuation occurs (Fig. 7a). Conversely, the result obtained from the new dynamic version shows consistency under different time steps (Fig. 7b). This clearly shows that the new dynamic version can use longer time steps while maintaining forecast accuracy, thereby greatly improving the integration efficiency of the model.

Figure 7.  The 850-hPa temperature (K) after integration over 9 days: (a) previous operational version and (b) new dynamic version. Black line: 600 s, red line: 1200 s, and green line: 1800 s.

Under horizontal resolution of 0.25°, the integration stability of both versions was tested for one month in winter and one month in summer. The previous operation version ensures stability with a 300-s time step, but approximately half of the case is unstable with a 450-s time step. The new dynamic version ensures stability for all cases with a 450-s time step, which represents an approximate one-third saving in integration time. The benefits realized in terms of computational stability are derived mainly from the use of the predictor–corrector scheme that avoids the computational noise caused by temporal extrapolation. The predictor–corrector scheme requires that the dynamic equations be solved twice per time step, but each process converges faster and by dramatically increasing the time step, the overall integration efficiency is increased by approximately 30%.

The CMA-GFS 3.0 has been operational in the CMA since June 2020. Before that, a 1-yr retrospective trial was conducted, and comprehensive evaluation in comparison with the previous operational version (CMA-GFS 2.4) was conducted. The retrospective trial from 1 December 2018 to 30 November 2019, covering four complete seasons, included four 6-h global 4DVAR analyses conducted daily and 10-day forecasts conducted at 0000 and 1200 UTC daily. The resolution was 0.25°, and the time step was 300 and 450 s for version 2.4 and version 3.0, respectively.

Through systematic testing and analysis, it can be seen that the overall prediction performance of version 3.0 is substantially enhanced in comparison with that of version 2.4. As shown in Fig. 8, the prediction capability regarding the circulation pattern has been systematically improved. The increase in the 500-hPa height anomaly correlation coefficient passed the 95% significance level for the first 6 days for both the Northern Hemisphere and the Southern Hemisphere, and it showed consistent improvement from the lower to upper levels.

Figure 8.  (a) Height anomaly correlation coefficient and (c) vertical profile for the Northern Hemisphere. (b, d) As in (a, c), respectively, but for Southern Hemisphere. Black line represents version 2.4; red line represents version 3.0.

As shown in Fig. 9, which compares the two groups of tests, the structure of the deviation in the temperature field is similar in the middle and lower layers, and changes obviously near 100 hPa; however, the performance of version 3.0 is better than that of version 2.4 in most regions. The deviations in the wind and humidity fields were found similar to those in the temperature field.

Figure 9.  Zonal mean temperature deviation (K) at 5 days: (a, c) version 2.4 and (b, d) version 3.0. (a, b) JJA and (c, d) DJF.

This paper introduces a new dynamic version of the CMA-GFS, which was developed in recent years, and compares the various aspects of performance between the new dynamic version and the previous operational version. The main conclusions derived are as follows.

(1) The core content of this version upgrade is the predictor–corrector SISL temporal integration scheme, which avoids temporal extrapolation of the upstream point and calculation of nonlinear terms. Furthermore, through combination with the 3D reference profile, the stability and accuracy are effectively improved, and the time discretization achieves second-order accuracy. At 0.25° resolution, the time step can be extended from 300 to 450 s, and the overall integration efficiency is increased by approximately 30%.

(2) Using the 3D reference profile, the reference state can be closer to the model atmosphere, and more information is included in the reference state without taking part in the SL discretization. This reduces the magnitude of the nonlinear terms and improves the accuracy of the spatial discretization. Currently, the 3D reference profile is set to the monthly average climate field, which satisfies the hydrostatic balance and does not change with time during the integration.

(3) With the hybrid vertical coordinate system, the coordinate plane of the upper level becomes flat, which can reduce the deviation of the PGF and scalar transport. With the new spectral filter terrain and optimization of artificial restrictions, improved distribution of the KES is obtained on the shortwave scale, which improves the simulation of small-scale circulations and terrain-forced precipitation.

(4) Benchmark tests, analysis of energy characteristics, and comprehensive evaluation of a 1-yr retrospective trial revealed that the overall prediction performance of the new dynamic version is substantially better than that of the previous operational version.

After more than 20 years of independent research and development, the CMA has established a complete operational numerical weather prediction system with the CMA-GFS model as the core, gradually narrowing the gap in terms of NWP capability to the ECMWF and NCEP. The new version dynamical core of the CMA-GFS has resolved the key problems of high-precision calculation and high-efficient integration in the autonomous model, providing a solid foundation for further development of the entire operational system of the CMA.

For future high-resolution applications on massively parallel computers, and for the foreseeable seamless prediction from weather to climate, substantial challenges remain regarding CMA-GFS dynamics. As with many models using implicit time-stepping, the nonlocality of the 3D elliptic equation solver becomes a critical issue regarding efficiency, and any nonlocal process will limit the speed of the model and the transposition of data. The anisotropy of the latitude–longitude grid will result in a larger number of solver iterations, particularly at high resolution. Therefore, more uniform grids on a spherical framework might be necessary as a replacement for the regular latitude–longitude grids.

Traditional SL schemes lose any conservation properties of the governing equations owing to point-wise interpolation of quantities to the upstream departure point. Although it is argued that mass conservation might not be that critical in weather forecasting, mass conservation should be revisited for future Earth system predictions, e.g., the seamless weather–climate modeling system and chemical weather forecasting. In the current CMA-GFS, conservation of moisture species is guaranteed by conservative SL scalar advection. Further development of the CMA-GFS model should consider dry air mass conservation by developing a conservative Lagrangian algorithm.

Physics coupling with 2TL SISL dynamics has become a hot topic among the NWP community. Dynamics–physics coupling is not carefully considered in the current CMA model. Introduction of a relatively optimal coupling algorithm would be very important in achieving a better balance between the physics forcing and the dynamics. This represents a critical issue for improving forecast skill.