
First, a brief introduction is given to the governing equations and the current temporal integration scheme of the CMAGFS. 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 CMAGFS 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:$$ \frac{{\rm{d}}u}{{\rm{d}}t}=\frac{{C}_{P}{\theta }_{v}}{r\mathit{{\rm{cos}}}\phi }\frac{\partial \Pi }{\partial \lambda }+fv , $$ (1) $$ \frac{{\rm{d}}v}{{\rm{d}}t}=\frac{{C}_{P}{\theta }_{v}}{r}\frac{\partial \Pi }{\partial \phi }fu , $$ (2) $$ {\delta }_{NH}\frac{{\rm{d}}w}{{\rm{d}}t}={C}_{P}{\theta }_{v}\frac{\partial \Pi }{\partial r}g , $$ (3) $$ (\gamma 1)\frac{{\rm{d }}\Pi }{{\rm{d}}t}=\Pi \cdot {D}_{3}+\frac{{F}_{\theta }^{*}}{{\theta }_{v}} , $$ (4) $$ \frac{{\rm{d}}\theta }{{\rm{d}}t}=\frac{{F}_{\theta }^{*}}{\Pi } , $$ (5) $$\begin{aligned}[b] & \frac{{\partial q}}{{\partial t}} + \frac{1}{{r\cos \varphi }}\frac{{\partial qu}}{{\partial \lambda }} + \frac{1}{r}\frac{{\partial qv}}{{\partial \varphi }} + \frac{{\partial qw}}{{\partial r}} = q\Bigg(\frac{1}{{r\cos \varphi }} \frac{{\partial u}}{{\partial \lambda }} \\ & + \frac{1}{r}\frac{{\partial v}}{{\partial \varphi }} + \frac{{\partial w}}{{\partial r}}\Bigg), \end{aligned} $$ (6) 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:
$$ \frac{{\rm{D}}\psi \left(X,t\right)}{{\rm{D}}t}=L\left(X,t,\psi \right)+N\left(X,t,\psi \right) , $$ (7) 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:$$ \frac{{\psi }^{n+1}{\psi }_{d}^{n}}{\Delta t}=\alpha {\left(L+N\right)}^{n+1}+\left(1\alpha \right){\left(L+N\right)}_{d}^{n} , $$ (8) where
$ n $ is the time level,$ \Delta t $ is the time step,$ d $ denotes the departure point, and$ \alpha $ is a timeweighted 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 midpoint$ \dfrac{{X}_{a}+{X}_{d}}{2} $ of the trajectory to obtain the wind fields:$$\hspace{20pt} {N}^{n+1}=2{N}^{n}{N}^{n1} , $$ (9) $$\hspace{20pt} V(X,{t}^{n+\frac{1}{2}})=\frac{3}{2}V\left(X,{t}^{n}\right)\frac{1}{2}V (X,{t}^{n1} ) . $$ (10) 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 firstorder 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 CMAGFS, an SI predictor–corrector scheme is proposed. As described in Section 2.2, a general Lagrangianform 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:
$$\begin{aligned}[b] & \frac{{\psi }^{*}{\psi }^{n}}{\Delta t}=\frac{3}{2}N\left({\psi }^{n}\right)\frac{1}{2}N\left({\psi }^{n1}\right),\\ & \frac{{\psi }^{n+1}{\psi }^{n}}{\Delta t}=\frac{1}{2}\left[N\left({\psi }^{*}\right)+N\left({\psi }^{n}\right)\right] , \end{aligned}$$ (11) which has secondorder 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 highresolution 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 CMAGFS, we use a modified secondorder Adams–Moulton formulation suggested by Durran and Blossey (2012) in the first step and trapezoidal averaging in the second step:
$$\begin{aligned}[b] & \frac{{\psi }^{*}{\psi }^{n}}{\Delta t}=\frac{3}{4}L\left({\psi }^{*}\right)+\frac{1}{4}L\left({\psi }^{n1}\right), \\ & \frac{{\psi }^{n+1}{\psi }^{n}}{\Delta t}=\frac{1}{2}\left[L\left({\psi }^{n+1}\right)+L\left({\psi }^{n}\right)\right]. \end{aligned} $$ (12) The nature of the stability of the above secondorder 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 twostep 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:
$$\begin{aligned} & \frac{{\psi }^{*}{{\psi }_{d}^{n}}}{\Delta t}= N{\left({\psi }^{*}\right)}^{n+\frac{1}{2}}+{\alpha }_{1}L\left({\psi }^{*}\right)+\left(1{\alpha }_{1}\right){L}_{d}\left({\psi }^{n1}\right), \\ & \frac{{\psi }^{n+1}{{\psi }_{d}^{n}}}{\Delta t}= \frac{1}{2}\left[N\left({\psi }^{*}\right)+{N}_{d}\left({\psi }^{n}\right)\right]+\frac{1}{2}\left[L\left({\psi }^{n+1}\right)+{L}_{d}\left({\psi }^{n}\right)\right] \\ & \quad\quad\quad = \; \left(1{\alpha }_{2}\right){\left(L+N\right)}_{d}^{n}+{\alpha }_{2}\left[N\left({\psi }^{*}\right)+L\left({\psi }^{n+1}\right)\right] ,\\[14pt] {} \end{aligned} $$ (13) where
$N{\left({\psi }^{*}\right)}^{n+\frac{1}{2}}=\dfrac{3}{2}N\left({\psi }^{n}\right)\dfrac{1}{2}N\left({\psi }^{n1}\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 CMAGFS, 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 CMAGFS, 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:
$$ {x}_{a}{x}_{d}\approx \Delta tV\left[x\left({t}^{n}+\frac{\Delta t}{2}\right),{t}^{n}+\frac{\Delta t}{2}\right] , $$ (14) 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:
$$ {x}_{d}={x}_{a}\frac{1}{2}\cdot \Delta t\cdot \left[{V}^{*}\left(\frac{{x}_{a}+{x}_{d}^{*}}{2}\right)+{V}^{n}\left(\frac{{x}_{a}+{x}_{d}}{2}\right)\right] , $$ (15) 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 CMAGFS 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 CMAGFS 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 righthand side of the discretized equations will include extra terms that take account of the horizontal variations of the reference state. 
Utilization of a terrainfollowing 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 terrainfollowing coordinate systems. In the first version of the CMAGFS, a heightbased terrainfollowing coordinate system was utilized but without flattening of the coordinate surfaces in the upper levels (GalChen and Somerville, 1975). Recently, a heightbased terrainfollowing hybrid vertical coordinate system was introduced to accompany implementation of the SI predictor–corrector scheme:
$$ z=a\Bigg({z}_{s}+\frac{{z}_{t}{z}_{s}}{{z}_{t}}\widehat{z}\Bigg)+(1a)\widehat{z},\; \left\{\begin{array}{l}z < {z}_{c1},\;a=1\\ {z}_{c1} < z < {z}_{c2}, \; 0 < a < 1\\ z\geqslant {z}_{c2},\;a=0\end{array}\right. , $$ (16) 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 CMAGFS 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:
$$\begin{aligned}[b] & {\psi }^{*}=\Delta t{\alpha }_{1}L\left({\psi }^{*}\right)+{{A}^{*}_{\psi }},\\ &{\rm here},\; {{A}^{*}_{\psi }}={{\psi }_{d}}^{ n}+\Delta t\Bigg[\left(1{\alpha }_{1}\right){L}_{d}\left({\psi }^{n1}\right)\\ & + \left(1+{\alpha }_{2}\right)N\left({\psi }^{n}\right){\alpha }_{2}N\left({\psi }^{n1}\right)\Bigg] .\end{aligned} $$ (17) 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 righthand 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 quasimonotonic SL scheme.
In the corrector step, the resulting scheme may be written as follows:
$$\begin{aligned}[b] & {\psi }^{n+1}=\Delta t{\alpha }_{2}L\left({\psi }^{n+1}\right)+{A}_{\psi}\\ & {\rm here},\hspace{0.33em}\hspace{0.33em} {A}_{\psi }={{\psi }_{d}^{n}}+\Delta t\left[{\alpha }_{2}N\left({\psi }^{*}\right)+\left(1{\alpha }_{2}\right){\left(L+N\right)}_{d}^{n}\right] .\end{aligned} $$ (18) Again, the only difference from the original scheme is in the second term of the righthand side of Eq. (18). It should be noted that the temporal weighting (or offcentered parameter)
$ {\alpha }_{2} $ takes a value of 0.5 in the corrector step to achieve secondorder 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 smallscale circulations and terrainforced precipitation. 2) With optimization of artificial restrictions, excessive damping and smoothing are removed, which improves simulation capabilities for smallscale 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 CMAGFS. The CMAGFS 3.0, which has been operational since June 2020, replaced the previous version of the model (i.e., CMAGFS 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.
CMAGFS 2.4
(previous operational version)CMAGFS 3.0
(new dynamic version)Temporal integration scheme Classical 2TL SISL Predictor–corrector 2TL SISL Reference profile 1D reference profile
based on isothermal atmosphere3D reference profile
based on climate fieldVertical coordinate system Heightbased terrainfollowing Hybrid heightbased terrainfollowing Terrain Relatively smooth terrain filter Higher precision terrain filter Artificial restriction 1) Fourthorder horizontal diffusion
2) Vertical velocity implicit damping (stronger)
3) Rayleigh friction at model top
4) Local humidity restrictions1) Vertical velocity implicit damping (weaker) Table 1. Comparison of CMAGFS 2.4 and CMAGFS 3.0
CMAGFS 2.4 (previous operational version)  CMAGFS 3.0 (new dynamic version)  
Temporal integration scheme  Classical 2TL SISL  Predictor–corrector 2TL SISL 
Reference profile  1D reference profile based on isothermal atmosphere  3D reference profile based on climate field 
Vertical coordinate system  Heightbased terrainfollowing  Hybrid heightbased terrainfollowing 
Terrain  Relatively smooth terrain filter  Higher precision terrain filter 
Artificial restriction  1) Fourthorder horizontal diffusion 2) Vertical velocity implicit damping (stronger) 3) Rayleigh friction at model top 4) Local humidity restrictions  1) Vertical velocity implicit damping (weaker) 