模型离散化

• Runge-Kutta 时间积分方案
• 声学积分
• Full Time-Split Integration Sequence
• Diabatic Forcing
• Hydrostatic Option

ARW 求解器使用时间分割积分方案。 一般而言，使用三阶 Runge-Kutta(RK3) 时间积分方案对慢速或低频（在气象上有意义的）模式进行积分，而在较小的时间步长上进行高频声学模式积分以保持数值稳定性。 使用前后时间积分方案对水平传播的声学模型（包括使用恒压上边界条件在质量坐标方程中存在的外部模型）和重力波进行积分，对垂直传播的声学模型与浮力振荡使用垂直隐式方案（使用声学时间步长）进行积分。

---

The time-split integration for the flux-form equations is described and analyzed in Klemp et al. (2007). The time-splitting is similar to that first developed by Klemp and Wilhelmson (1978) for leapfrog time integration and analyzed by Skamarock and Klemp (1992). This time-split approach was extended to the RK3 scheme as described in Wicker and Skamarock (2002).

---

参考文献中描述的较早实现方式与 ARW 实现方式之间的主要区别与我们使用质量垂直坐标和通量形式的方程组有关，如 Klemp 等人 (2007 年）所述，以及我们在时间分割积分的声学成分中使用摄动变量。 声学模型积分以对 RK3 积分的校正形式进行转换。

Runge-Kutta 时间积分方案

Wicker 和 Skamarock(2002) 中描述的 RK3 方案使用 predictor-corrector 公式集成了一组常微分方程。 定义 ARW 求解器中的预报变量为 $\Phi = (U,V,W,\Theta_m,\phi',\mu'_d,Q_m)$，定义模型方程为 $\Phi_t = R(\Phi)$，RK3 积分采取 3 步形式来推进，从 $\Phi(t)$ 到 $\Phi(t + \Delta t)$ 的解：

$$\Phi^* = \Phi^t + \frac{\Delta t}{3} R(\Phi^t) \tag{3.1}$$

$$\Phi^{**} = \Phi^t + \frac{\Delta t}{2} R(\Phi^*) \tag{3.2}$$

$$\Phi^{t + \Delta t} = \Phi^t + \Delta t R(\Phi^{**}) \tag{3.3}$$

其中 $\Delta t$ 是低频模式的时间步长（模型时间步长）。 在 (3.1) – (3.3) 中，上标表示时间级别。 该方案本身并不是真正的 Runge-Kutta 方案，因为尽管它对于线性方程是三阶精确的，但对于非线性方程只是二阶的。 对于 ARW 方程，时间导数 $\Phi_t$ 是方程 (2.21)，(2.24) 和 (2.28) – (2.32) 中的偏时间导数（最左边的项），而 $R(\Phi)$ 是方程中的其余项。

声学积分

在气象意义上微不足道的高频声学模式将严重限制 (3.1) – (3.3) 中的 RK3 时间步长 $\Delta t$。为了避免这种时间步长的限制，我们使用了 Wicker 和 Skamarock(2002) 中描述的时间分割方法。 另外，为了提高分割的准确性，我们在 RK3 长时间步长序列中使用较小的声学时间步长来合并控制方程的扰动形式。 为了塑造 RK3 时间分割声学积分的扰动方程，我们定义了较小的时间步长变量，这些变量与最新的 RK3 预测变量（由上标 $t^*$ 表示，也表示 (3.1) - (3.3) 中的 $\Phi^t, \Phi^∗$ 或 $\Phi^{**}$）存在偏差：

$$\begin{aligned} V'' &= V - V^{t^*},& \Omega'' &= \Omega - \Omega^{t^*},& \Theta_m'' &= \Theta_m - \Theta_m^{t^*} \\ \phi'' &= \phi' - {\phi'}^{t^*},& \alpha_d'' &= \alpha_d' - {\alpha_d'}^{t^*},& \mu_d'' &= \mu_d' - {\mu_d'}^{t^*} \end{aligned}$$

静力关系（即垂直坐标定义）变为

$$\alpha''_d = - \frac{1}{\mu^{t*}_d} (\partial_{\eta}\phi'' + \alpha''_d \mu''_d) \tag{3.4}$$

此外，我们还介绍了状态方程的一个版本，该版本对 $t^*$ 线性化，

$$p'' = \frac{c_s^2}{\alpha_d^{t*}} (\frac{\Theta_m''}{\Theta_m^{t*}} - \frac{\alpha_d''}{\alpha_d^{t*}} - \frac{\mu_d''}{\mu_d^{t*}}) \tag{3.5}$$

其中 $c^2_s=\gamma p^{t*} \alpha^{t*}_d$ 是声速的平方。 根据模型的预测变量，使用线性状态方程 (3.5) 和垂直坐标定义 (3.4) 来转换 (2.30) 中的垂直压力梯度。 通过组合 (3.5) 和 (3.4)，垂直压力梯度可以表示为

$$\partial_\eta p'' = \partial_\eta (C \partial_\eta \phi'') +\partial_\eta (\frac{c_s^2}{\alpha_d^{t*}} \frac{\Theta_m''}{\Theta_m^{t*}}) \tag{3.6}$$

其中 $C=c^2_s / {\mu_d^{t*}}^2$ 。 关于当前大时间步长的线性化应该在几个小时间步长的时间间隔内高度准确。

这些变量与 (3.6) 一起代入方程 (2.21) 和 (2.28) – (2.33) ，并得出声学时步方程：

$$\begin{aligned} &\partial_{t} U^{\prime \prime} + \\ &\left(m_{x} / m_{y}\right)\left(\alpha^{t^{*}} / \alpha_{d}^{t^{*}}\right)\left[\mu_{d}^{t^{*}}\left(\alpha_{d}^{t^{*}} \partial_{x} p^{\prime \prime \tau}+\alpha_{d}^{\prime \prime \tau} \partial_{x} \bar{p}+\partial_{x} \phi^{\prime \prime \tau}\right)+\partial_{x} \phi^{t^*}\left(\partial_{\eta} p^{\prime \prime}-\mu_{d}^{\prime \prime}\right)^{\tau}\right]=R_{U}^{t^{*}} \tag{3.7} \end{aligned}$$

$$\begin{aligned} &\partial_{t} V^{\prime \prime} + \\ &\left(m_{y} / m_{x}\right)\left(\alpha^{t^{*}} / \alpha_{d}^{t^{*}}\right)\left[\mu_{d}^{t^{*}}\left(\alpha_{d}^{t^{*}} \partial_{y} p^{\prime \prime \tau}+\alpha_{d}^{\prime \prime \tau} \partial_{y} \bar{p}+\partial_{y} \phi^{\prime \prime \tau}\right)+\partial_{y} \phi^{t^*}\left(\partial_{\eta} p^{\prime \prime}-\mu_{d}^{\prime \prime}\right)^{\tau}\right]=R_{V}^{t^{*}} \tag{3.8} \end{aligned}$$

$$\delta_{\tau} \mu_{d}^{\prime \prime}+m_{x} m_{y}\left[\partial_{x} U^{\prime \prime}+\partial_{y} V^{\prime \prime}\right]^{\tau+\Delta \tau}+m_{y} \partial_{\eta} \Omega^{\prime \prime \tau+\Delta \tau}=R_{\mu}^{t^{*}} \tag{3.9}$$

$$\delta_{\tau} \Theta_{m}^{\prime \prime}+m_{x} m_{y}\left[\partial_{x}\left(U^{\prime \prime} \theta_{m}^{t^{*}}\right)+\partial_{y}\left(V^{\prime \prime} \theta_{m}^{t^{*}}\right)\right]^{\tau+\Delta \tau}+m_{y} \partial_{\eta}\left(\Omega^{\prime \prime \tau+\Delta \tau} \theta_{m}^{t^*}\right)=R_{\Theta_{m}}^{t^{*}} \tag{3.10}$$

$$\delta_{\tau} W^{\prime \prime}-m_{y}^{-1} \overline{g\left\{\left(\alpha / \alpha_{d}\right)^{t^{*}}\left[\partial_{\eta}\left(C \partial_{\eta} \phi^{\prime \prime}\right)+\partial_{\eta}\left(\frac{c_{s}^{2}}{\alpha^{t}} \frac{\Theta_{m}^{\prime \prime}}{\Theta_{m}^{t^{\prime}}}\right)\right]-\mu_{d}^{\prime \prime}\right\}}^{\prime}=R_{W}^{t^{*}} \tag{3.11}$$

$$-\delta_{\tau} \phi^{\prime \prime}+\frac{1}{\mu_{d}^{t^{*}}}\left[m_{y} \Omega^{\prime \prime \tau+\Delta \tau} \delta_{\eta} \phi^{t^{*}}-m_{y} \overline{g W^{\prime \prime}}^{\prime \prime}\right]=R_{\phi}^{t^{*}} \tag{3.12}$$

(3.7) – (3.12) 中的右侧项对于包含每个 RK3 子步骤（即 (3.1) – (3.3)) 的时间积分的声学步骤是固定的，并且由

$$\begin{aligned} R_{U}^{t^{*}}=&-m_{x}\left[\partial_{x}(U u)+\partial_{y}(V u)\right]-\partial_{\eta}(\Omega u) \\ &-\left(m_{x} / m_{y}\right)\left(\alpha / \alpha_{d}\right)\left[\mu_{d}\left(\partial_{x} \phi^{\prime}+\alpha_{d} \partial_{x} p^{\prime}+\alpha_{d}^{\prime} \partial_{x} \bar{p}\right)+\partial_{x} \phi\left(\partial_{\eta} p^{\prime}-\mu_{d}^{\prime}\right)\right] \tag{3.13} \end{aligned}$$

$$\begin{aligned} R_{V}^{t^{*}}=&-m_{y}\left[\partial_{x}(U v)+\partial_{y}(V v)\right]-\left(m_{y} / m_{x}\right) \partial_{\eta}(\Omega v) \\ &-\left(m_{y} / m_{x}\right)\left(\alpha / \alpha_{d}\right)\left[\mu_{d}\left(\partial_{y} \phi^{\prime}+\alpha_{d} \partial_{y} p^{\prime}+\alpha_{d}^{\prime} \partial_{y} \bar{p}\right)+\partial_{y} \phi\left(\partial_{\eta} p^{\prime}-\mu_{d}^{\prime}\right)\right] \tag{3.14} \end{aligned}$$

$$R_{\mu_{d}}^{t^{*}}=-m_{x} m_{y}\left[\partial_{x} U+\partial_{y} V\right]-m_{y} \partial_{\eta} \Omega \tag{3.15}$$

$$R_{\Theta_{m}}^{t^{*}}=-m_{x} m_{y}\left[\partial_{x}\left(U \theta_{m}\right)+\partial_{y}\left(V \theta_{m}\right)\right]-m_{y} \partial_{\eta}\left(\Omega \theta_{m}\right)+F_{\Theta_{m}} \tag{3.16}$$

$$\begin{aligned} R_{W}^{t^{*}}=&-m_{x}\left[\partial_{x}(U w)+\partial_{y}(V w)\right]-\partial_{\eta}(\Omega w) \\ &+m_{y}^{-1} g\left(\alpha / \alpha_{d}\right)\left[\partial_{\eta} p^{\prime}-\bar{\mu}_{d}\left(q_{v}+q_{c}+q_{r}\right)\right]-m_{y}^{-1} \mu_{d}^{\prime} g+F_{W} \tag{3.17} \end{aligned}$$

$$R_{\phi}^{t^{*}}=-\mu_{d}^{-1}\left[m_{x} m_{y}\left(U \partial_{x} \phi+V \partial_{y} \phi\right)+m_{y} \Omega \partial_{\eta} \phi-m_{y} g W\right] \tag{3.18}$$

其中 (3.13) – (3.18) 中的所有变量均在时间 $t^*$ 进行评估（即，对于 (3.1) – (3.3) 中的 RK3 子步骤，使用 $\Phi^t，\Phi^* 或 \Phi^{**}$ )。公式 (3.7) – (3.12) 利用离散声学时间步长算子

$$\delta_{\tau} a=\frac{a^{\tau+\Delta \tau}-a^{\tau}}{\Delta \tau}$$

其中 $\Delta \tau$ 是声学时间步长，同时平均算子将在声学时间步长的时间平均项上稍微向前倾斜

$$\bar{a}^{\tau}=\frac{1+\beta}{2} a^{\tau+\Delta \tau}+\frac{1-\beta}{2} a^{\tau} \tag{3.19}$$

其中，$\beta$ 是用户指定的参数（请参见第 4.3.3 节）。

声学时间步长上的积分如下进行。 从时间 $\tau$ 处以小时间步长变量开始，向前推（3.7）和（3.8）以获得 $U''^{\tau+\Delta\tau}$ 和 $V''^{\tau+\Delta\tau}$。 然后从（3.9）计算出 $\mu''^{\tau+\Delta\tau}$ 和 $\Omega''^{\tau+\Delta\tau}$。 这是通过首先从地表到区域顶部垂直积分（3.9）来完成的，这消除了 $\partial_\eta \Omega''$。 回顾 $\mu_d = \partial p_d / \partial \eta$ 且混合垂直坐标的 $p_d$ 由（2.2）定义，则（3.9）的垂直积分变为，

$$\delta_{\tau} p_{c}^{\prime \prime}=m_{x} m_{y} \int_{1}^{0}\left[\partial_{x} U^{\prime \prime}+\partial_{y} V^{\prime \prime}\right]^{\tau+\Delta \tau} d \eta \tag{3.20}$$

其中，$p_c(x，y) = p_s-p_t$ 是（x，y）处垂直气柱中的干静力压力差（质量）。 从（3.20）计算出 $\delta_\tau p''_c$ 之后，通过使用地表的下边界条件 $\Omega'' = 0$ 对（3.9）项 $\partial_\eta \Omega''$ 进行垂直积分（$\delta_\tau \mu = B_\eta \delta_\tau p_c$）来获得 $\Omega''^{\tau+\Delta\tau}$，得

$$\Omega^{\prime \prime \tau+\Delta \tau}=m_{y}^{-1}[1-B(\eta)] \delta_{\tau} p_{c}-m_{x} \int_{1}^{\eta}\left[\partial_{x} U^{\prime \prime}+\partial_{y} V^{\prime \prime}\right]^{\tau+\Delta \tau} d \eta \tag{3.21}$$

and µ 00τ+∆τ d is recovered using (2.2):

$$\mu_{d}^{\prime \prime \tau+\Delta \tau}(x, y, \eta)=B_{\eta}(\eta) p_{c}^{\prime \prime \tau+\Delta \tau}(x, y)+\left[1-B_{\eta}(\eta)\right]\left(p_{0}-p_{t}\right) \tag{3.22}$$

From (3.22), it is evident that µ need not be stored as a three-dimensional array, but can be readily constructed when needed from the two-dimensional p c (x,y) array together with the one-dimensional B η (η) profile.

Knowing Ω 00 τ+∆τ , (3.10) can be stepped forward to calculate Θ 00 m τ+∆τ . Equations (3.11) and (3.12) are combined to form a vertically implicit equation that is solved for W 00 τ+∆τ subject to the boundary condition Ω = Ω 00 = 0 at the surface (z = h(x,y)) and p 0 = 0 along the model top. φ 00 τ+∆τ is then obtained from (3.12), and p 00 τ+∆τ and α 00 d τ+∆τ are recovered from (3.5) and (3.4).

Full Time-Split Integration Sequence

The time-split RK3 integration technique is summarized in Figure 3.1. It consists of two primary loops— an outer loop for the large-time-step Runge-Kutta integration, and an inner loop for the acoustic mode integration.

In the RK3 scheme, physics can be integrated within the RK3 time integration (using a time forward step, i.e., step (1) in Fig. 3.1, or the RK3 time integration if higher temporal accuracy is desired, i.e., in step (2)— implying a physics evaluation every RK3 substep) or external to it using additive timesplitting, i.e., step (9).

Within the acoustic integration, the acoustic time step ∆τ is specified by the user through the choice of n s (see Section 3.3.2). Within the first RK3 substep, however, a single acoustic time step is used to advance the solution regardless of n s . Within the full RK3-acoustic timesplit integration, this modified acoustic time step does not impose any additional stability constraints (see Wicker and Skamarock, 2002).

The major costs in the model arise from the evaluation of the right hand side terms R t ∗ in (3.7) – (3.12). The efficiency of the RK3 timesplit scheme arises from the fact that the RK3 time step ∆t is much larger than the acoustic time step ∆τ, hence the most costly evaluations are only performed in the less-frequent RK3 steps.

Diabatic Forcing

Within the RK3 integration sequence outlined in Fig. 3.1, the RHS term R t ∗ Θ m in the thermo- dynamic equation (3.10) contains contributions from the diabatic physics tendencies that are computed in step (1) at the beginning of the first RK3 step. This diabatic forcing is integrated within the acoustic steps (specifically, in step 4 in the time integration sequence shown in Fig.

Begin Time Step
    Begin RK3 Loop: Steps 1, 2, and 3
    (1) If RK3 step 1, compute and store F Φ
        (i.e., physics tendencies for RK3 step, including mixing).
    (2) Compute R t ∗ , (3.13)–(3.18)
    Begin Acoustic Step Loop: Steps 1 → n,
        RK3 step 1, n = 1, ∆τ = ∆t/3;
        RK3 step 2, n = n s /2, ∆τ = ∆t/n s ;
        RK3 step 3, n = n s , ∆τ = ∆t/n s .
        (3) Advance horizontal momentum, (3.7) and (3.8)
            Global: Apply polar filter to U 00τ+∆τ and V 00τ+∆τ .
        (4) Advance µ d (3.9) and compute Ω 00 τ+∆τ then advance Θ m (3.10)
            Global: Apply polar filter to µ τ+∆τ d and Θ 00τ+∆τ m .
        (5) Advance W and φ (3.11) and (3.12)
            Global: Apply polar filter to W 00τ+∆τ and φ 00τ+∆τ .
        (6) Diagnose p 00 and α 00 using (3.5) and (3.4)
    End Acoustic Step Loop
    (7) Scalar transport: Advance scalars (2.24)
        over RK3 substep (3.1), (3.2) or (3.3)
        (using mass fluxes U, V and Ω time-averaged over the acoustic steps).
        Global: Apply polar filter to scalars.
    (8) Using updated prognostic variables, compute p 0 with (2.16) and α 0 with (2.33) or (2.34)
    End RK3 Loop
    (9) Compute non-RK3 physics (currently microphysics), advance variables.
    Global: Apply polar filter to updated variables.
End Time Step

Figure 3.1: Time step integration sequence. Here n represents the number of acoustic time steps for a given substep of the RK3 integration, and n s is the ratio of the RK3 time step to the acoustic time step for the second and third RK3 substeps.

3.1). Additional diabatic contributions are integrated in an additive-time-split manner in step (9) after the RK3 update is complete. Thus, the diabatic forcing computed in step (9) (the microphysics in the current release of the ARW) does not appear in R t ∗ Θ m from (3.10) in the acoustic integration. We have discovered that this time splitting can excite acoustic waves and can give rise to noise in the solutions for some applications. Note that the non-RK3 physics are integrated in step (9) because balances produced in the physics are required at the end of the time step (e.g., the saturation adjustment in the microphysics). So while moving these non-RK3 physics into step (1) would eliminate the noise, the balances produced by these p hysics would be altered.

We have found that the excitation of the acoustic modes can be circumvented while leaving the non-RK3 physics in step (9) by using the following procedure that is implemented in the ARW. In step (1) of the integration procedure (Fig. 3.1), an estimate of the diabatic forcing in the Θ m equation arising from the non-RK3 physics in step (9) is included in the diabatic forcing term R t ∗ Θ m in (3.10) (which is advanced in step 4). This estimated diabatic forcing is then removed from the updated Θ m after the RK3 integration is complete and before the evaluation of the non-RK3 physics in step (9). We use the diabatic forcing from the previous time step as the estimated forcing; hence this procedure results in few additional computations outside of saving the diabatic forcing between time steps.

Hydrostatic Option

A hydrostatic option is available in the ARW solver. The time-split RK3 integration technique summarized in Fig. 3.1 is retained, including the acoustic step loop. Steps (5) and (6) in the acoustic-step loop, where W and φ are advanced and p 00 and α 00 are diagnosed, are replaced by the following three steps. (1) Diagnose the pressure from the full hydrostatic equation (including moisture)

$$\delta_{\eta} p_{h}=\frac{\alpha_{d}}{\alpha} \mu_{d}=\left(1+\sum q_{m}\right) \mu_{d}$$

(2) Diagnose α d using the equation of state (2.16) and the prognosed θ m . (3) Diagnose the geopotential using the dry hydrostatic equation (2.33). The vertical velocity w can be diagnosed from the geopotential equation, but it is not needed in the solution procedure. The acoustic step loop advances gravity waves, including the external mode, and the Lamb wave when the hydrostatic option is used.