﻿ 全球质量守恒准均匀经纬网格三次样条函数变换准拉格朗日积分方案与模拟个例
 文章快速检索 高级检索
 高原气象  2017, Vol. 36 Issue (4): 1091-1105  DOI: 10.7522/j.issn.1000-0534.2016.00069 0

### 引用本文 [复制中英文]

[复制中文]
Gu Xuzan, Zhao Jun, Tang Yonglan. 2017. A Quasi-Lagrangian Integration of Conservation of Atmospheric Mass with Unify Scheme of Cubic Spline Function Transformating on Quasi-uniform Latitude-longitude Grid and Its Integration Cases[J]. Plateau Meteorology, 36(4): 1091-1105. DOI: 10.7522/j.issn.1000-0534.2016.00069.
[复制英文]

### 文章历史

1. 中国气象局武汉暴雨研究所暴雨监测预警湖北省重点实验室, 武汉 430074;
2. 中国人民解放军国防科学技术大学海洋科学与工程研究院, 长沙 410073

1 引言

2 样条格式准拉格朗日预报方程 2.1 大气运动方程

 $\frac{{{\rm{d}}u}}{{{\rm{d}}t}} = - RT\frac{{\partial \ln p}}{{\partial x}} + fv - \tilde fw + \frac{{uv \cdot \tan \varphi - uw}}{{{r_0}}} + {F_u}\hat = {a_u},$ (1)
 $\frac{{{\rm{d}}v}}{{{\rm{d}}t}} = - RT\frac{{\partial \ln p}}{{\partial y}} - fu - \frac{{{u^2}\tan \varphi + vw}}{{{r_0}}} + {F_v}\hat = {a_v},$ (2)
 $\frac{{{\rm{d}}w}}{{{\rm{d}}t}} = - RT\frac{{\partial \ln p}}{{\partial z}} - g + \tilde fu + \frac{{{u^2} + {v^2}}}{{{r_0}}} + {F_w}\hat = {a_w},$ (3)
 $\frac{{\partial \ln p}}{{\partial t}} = \frac{{ - 1}}{{1 - \kappa }}\nabla \cdot V\hat = {a_p},$ (4)
 $\frac{{\partial \ln T}}{{\partial t}} = \frac{{ - \kappa }}{{1 - \kappa }}\nabla \cdot V\hat = {a_T},$ (5)
 $\frac{{{\rm{d}}q}}{{{\rm{d}}t}} = 0,$ (6)

uu方程在极点无定义, 而vv方程在极点有定义。先定义北、南极点(分别用下标NS表示)与0°E u平行(λ=0) 分量为uNuS, 而与0°E v平行分量(λ=π9/2) 为vNvS, 则对于(uN, vN)和(uS, vS), (任何水平矢量)都可以作如下三角函数矢量分解:

 $\left\{ \begin{array}{l} {u_{N\left( \lambda \right)}} = {u_N}\cos \lambda + {v_N}\sin \lambda \\ {u_{S\left( \lambda \right)}} = {u_S}\cos \lambda - {v_S}\sin \lambda \\ {v_{N\left( \lambda \right)}} = {v_N}\cos \lambda - {u_N}\sin \lambda \\ {v_{S\left( \lambda \right)}} = {v_S}\cos \lambda - {u_S}\sin \lambda \end{array} \right.,$ (7)

 $\frac{{{\rm{d}}{u_N}}}{{{\rm{d}}t}} = - RT{\left( {\frac{{\partial \ln p}}{{\partial y}}} \right)_{N\left( {3{\rm{\pi /2}}} \right)}} + f{v_N} - \frac{{{u_N}w}}{{{r_0}}} + {F_v},$ (8)
 $\frac{{{\rm{d}}{u_S}}}{{{\rm{d}}t}} = - RT{\left( {\frac{{\partial \ln p}}{{\partial y}}} \right)_{S\left( {{\rm{\pi /2}}} \right)}} + f{v_S} - \frac{{{u_S}w}}{{{r_0}}} + {F_v},$ (9)
 $\frac{{{\rm{d}}{v_N}}}{{{\rm{d}}t}} = - RT{\left( {\frac{{\partial \ln p}}{{\partial y}}} \right)_{N\left( 0 \right)}} + f{u_N} - \frac{{{v_N}w}}{{{r_0}}} + {F_v},$ (10)
 $\frac{{{\rm{d}}{v_S}}}{{{\rm{d}}t}} = - RT{\left( {\frac{{\partial \ln p}}{{\partial y}}} \right)_{S\left( {\rm{0}} \right)}} - f{u_S} - \frac{{{v_S}w}}{{{r_0}}} + {F_v}.$ (11)
2.2 准拉格朗日预报方程

P表示pTquvw等, 已知P的一阶时间全导数(变率)为$\frac{{{\rm{d}}P}}{{{\rm{d}}t}}\hat = a$, 对于uvw方程, (au, av, aw)为旋转地球上单位质量空气所受广义牛顿力(加速度), 且它们是“气压梯度力+重力+地转偏向力+曲率力+摩擦力”之合力的三个分力; 对于pT方程, (ap, aT)为三维散度造成的绝热变率; 对于q方程, 其为水汽源汇与物理相变变率, 干绝热时其为零。

 $\begin{array}{l} P\left( {t + \Delta t,x,y,z} \right) = P\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + a\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\Delta t\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + a'\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\frac{{\Delta {t^2}}}{2}, \end{array}$ (12)

 $\begin{array}{l} P\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right) \approx P\left( {t,x,y,z} \right) - \Delta x\frac{{\partial P}}{{\partial x}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; - \Delta y\frac{{\partial P}}{{\partial y}} - \Delta z\frac{{\partial P}}{{\partial z}} + \frac{{{\Delta ^2}x}}{2}\frac{{{\partial ^2}P}}{{\partial {x^2}}} + \frac{{{\Delta ^2}y}}{2}\frac{{{\partial ^2}P}}{{\partial {y^2}}} + \frac{{{\Delta ^2}z}}{2}\frac{{{\partial ^2}P}}{{\partial {z^2}}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; + \Delta x\Delta y\frac{{{\partial ^2}P}}{{\partial x\partial y}} + \Delta x\Delta z\frac{{{\partial ^2}P}}{{\partial x\partial z}} + \Delta y\Delta z\frac{{{\partial ^2}P}}{{\partial y\partial z}}, \end{array}$ (13)

 $\left\{ \begin{array}{l} \Delta x = u\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\Delta t\\ \;\;\;\;\;\;\;\; + {a_u}\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\frac{{{\Delta ^2}t}}{2}\\ \Delta y = v\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\Delta t\\ \;\;\;\;\;\;\;\; + {a_v}\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\frac{{{\Delta ^2}t}}{2}\\ \Delta z = w\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\Delta t\\ \;\;\;\;\;\;\;\; + {a_w}\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\frac{{{\Delta ^2}t}}{2} \end{array} \right.,$ (14)

2.3 样条格式二阶时空离散预报方程

 $\begin{gathered} P\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right) \approx P\left( {t,x,y,z} \right) - \Delta x{P^x} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; - \Delta y{P^y} - \Delta z{P^z} + \frac{{\Delta {x^2}}}{2}{P^{xx}} + \frac{{\Delta {y^2}}}{2}{P^{yy}} + \frac{{\Delta {z^2}}}{2}{P^{zz}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; + \Delta x\Delta y{P^{xy}} + \Delta x\Delta z{P^{xz}} + \Delta y\Delta z{P^{yz}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; = \dddot P\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right), \hfill \\ \end{gathered}$ (15)

 $\begin{gathered} \dddot P\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right) \approx P\left( {t,x,y,z} \right) - \Delta x{P^x} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; - \Delta y{P^y} - \Delta z{P^z} + \frac{{\Delta {x^2}}}{2}{P^{xx}} + \frac{{\Delta {y^2}}}{2}{P^{yy}} + \frac{{\Delta {z^2}}}{2}{P^{zz}} + \Delta x\Delta y{P^{xy}} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; = \ddot P\left( {t,x - \Delta x,y - \Delta y,z} \right) - \Delta z{P^z} + \frac{{{\Delta ^2}z}}{2}{P^{zz}}, \hfill \\ \end{gathered}$ (16)

 $\begin{gathered} P\left( {t + \Delta t,x,y,z} \right) = \dddot P\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right) \hfill \\ \;\;\;\;\;\;\;\;\; + \dddot a\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\Delta t \hfill \\ \;\;\;\;\;\;\;\;\; + \dddot a'\left( {t,x - \Delta x,y - \Delta y,z - \Delta z} \right)\frac{{\Delta {t^2}}}{2}, \hfill \\ \end{gathered}$ (17)

 ${P^{t + \Delta t}} = \dddot P + \dddot a\Delta t + \dddot a'\frac{{\Delta {t^2}}}{2},$ (18)

 ${P^{t + \Delta t}} = \dddot P + \dddot a\Delta t,$ (19)

 ${P^{t + \Delta t}} = \dddot P + \dddot a\Delta t + c\frac{{\Delta {t^2}}}{2},$ (20)

 ${a^{t + \Delta t}} = \dddot a + \dddot a'\Delta t = \dddot ac\Delta t,$ (21)

 ${P^{t + \Delta t}} = \dddot P + \left( {{a^{t + \Delta t}} + \dddot a} \right)\frac{{\Delta t}}{2},$ (22)

3 样条格式与中央差比较 3.1 非线性样条格式

 $\begin{array}{l} {{\bar m}_j} = {m_j} + \frac{{{m_{j + 1}} - 2{m_j} + {m_{j - 1}}}}{6}\\ \;\;\;\;\; = \frac{{{P_{j + 1}} - {P_{j - 1}}}}{{2\Delta X}}, \end{array}$ (23)
 $\begin{array}{l} {{\bar M}_j} = {M_j} + \frac{{{M_{j + 1}} - 2{M_j} + {M_{j - 1}}}}{6}\\ \;\;\;\;\; = \frac{{{P_{j + 1}} - 2{P_j} + {P_{j - 1}}}}{{\Delta {X^2}}}, \end{array}$ (24)

3.2 样条格式线性截断误差

 $\frac{{\partial P}}{{\partial x}} = ik\exp \left( {ikx} \right),$ (25)
 $\frac{{{\partial ^2}P}}{{\partial {x^2}}} = - {k^2}\exp \left( {ikx} \right),$ (26)

x=jΔX, 可以用一阶中央差分mj近似表示一阶偏微商:

 $\begin{array}{l} {{\bar m}_j} = \frac{{{P_{j + 1}} - {P_{j - 1}}}}{{2\Delta X}}\\ \;\;\;\;\; = \frac{{\exp \left[ {ik\left( {j + 1} \right)\Delta X} \right] - \exp \left[ {ik\left( {j - 1} \right)\Delta X} \right]}}{{2\Delta X}}\\ \;\;\;\;\; = \frac{{i\sin k\Delta X}}{{\Delta X}}\exp \left( {ikj\Delta X} \right), \end{array}$ (27)

 ${k_m} = \frac{{\sin k\Delta X}}{{\Delta X}},$ (28)

 $\begin{array}{l} {{\bar M}_j} = \frac{{{P_{j + 1}} - 2{P_j} + {P_{j - 1}}}}{{\Delta {X^2}}}\\ \;\;\;\;\; = \frac{{\exp \left[ {ik\left( {j + 1} \right)\Delta X} \right] - 2\exp \left( {ikj} \right)\Delta X + \exp \left[ {ik\left( {j - 1} \right)\Delta X} \right]}}{{\Delta {X^2}}}\\ \;\;\;\;\; = - \frac{{2\left( {1 - \cos k\Delta X} \right)}}{{\Delta {X^2}}}\exp \left( {ikj\Delta X} \right), \end{array}$ (29)

 ${k_M} = \frac{{{{\left[ {2\left( {1 - \cos k\Delta X} \right)} \right]}^{\frac{1}{2}}}}}{{\Delta X}},$ (30)

${R_m} = \frac{{{k_m}}}{k} = \frac{{{\rm{sin}}k\Delta X}}{{k\Delta X}} \le 1,{R_M} = \frac{{{k_M}}}{k} = \frac{{{{\left[ {2\left( {1 - {\rm{cos}}k\Delta X} \right)} \right]}^{\frac{1}{2}}}}}{{k\Delta X}} \le 1$, RmRM可以反映它们的空间截断误差(比值):当RmRM接近于1时, 其能够比较准确地拟合简谐波及其空间微商(表 1)。

3.3 样条格式线性相速、群速与误差

 $\frac{{\partial P}}{{\partial t}} + {u_0}\frac{{\partial P}}{{\partial x}} = 0,$ (31)

$\frac{\partial }{{\partial t}}$式(31) 运算, 有:

 $\frac{{{\partial ^2}P}}{{\partial {t^2}}} - u_0^2\frac{{{\partial ^2}P}}{{\partial {x^2}}} = 0,$ (32)

 $i\left( { - \omega + {u_0}\frac{{\sin k\Delta X}}{{\Delta X}}} \right)\exp \left[ {i\left( {kj\Delta X - \omega t} \right)} \right] = 0,$ (33)
 $\left[ { - {\omega ^2} + u_0^2\frac{{2\left( {1 - \cos k\Delta X} \right)}}{{\Delta {X^2}}}} \right]\exp \left[ {i\left( {kj\Delta X - \omega t} \right)} \right] = 0,$ (34)

 ${c_m} = \frac{\omega }{k} = {u_0}\frac{{\sin k\Delta X}}{{k\Delta X}},$ (35)
 ${C_m} = \frac{{\partial \omega }}{{\partial k}} = {u_0}\cos k\Delta X,$ (36)
 ${c_M} = {u_0}\frac{{{{\left[ {2\left( {1 - \cos k\Delta X} \right)} \right]}^{\frac{1}{2}}}}}{{k\Delta X}},$ (37)
 ${C_M} = {u_0}\frac{{\sqrt 2 \sin k\Delta X}}{{2{{\left( {1 - \cos k\Delta X} \right)}^{\frac{1}{2}}}}},$ (38)

4 准均匀经纬网格设计 4.1 球面经纬网格

1°×1° A网格共有(0:360, 0:180)65160个网格点。为采用周期三次样条做拟合, 在纬向取0°E/0°W为重合点, 即设定变量场在起始点‘0’(0°E)值始终等于终止点‘360’(0°W)值; 并且标量场(如pTq)及垂直矢量场(如w)在两极点始终赋予360个相同值, 而水平矢量场(如(u, v)、(au, av))始终赋予360个不同“三角函数”分解值(参见式(7))。

4.2 准均匀经纬网格

(1) A-B流程(A-B网格切换): ① 在A网格上, 做变量场双三次曲面拟合(辜旭赞, 2011), 实现球面标、矢量场二阶可导, 以对B网格点(包含极点)做水平平流预报。② 在B网格上, 做水平平流变量和位移纬向三次样条拟合, 从而给全部A网格点赋予“预报”(水平平流和位移)插值。③ 在A网格上, 做“水平平流位移”经向和纬向三次样条拟合, 求出样条格式水平散度场, 以对B网格点做平流增压和绝热增温预报, (并可做气块水汽相变及降水预报)。④ 在B网格上, 做全部水平平流变量垂直三次样条拟合, 从而对B网格点做静力垂直平流预报。⑤ 在B网格上, 再做全部“水平平流+垂直平流”变量纬向三次样条拟合, 以给全部A网格点赋予“预报”插值。

(2) A-J流程(A-J网格切换): ①-⑤ 步骤均同于“A-B流程”, 但将其中的“B网格”通换成“J网格”。

(3) A(B)-J流程(A(B)-J网格切换): ① 在A网格上, 做变量场双三次曲面拟合, 实现球面标、矢量场二阶可导, 以对J网格点(包含极点)做水平平流预报。② 在J网格上, 做水平平流变量和位移纬向三次样条拟合, 以给全部A网格点(包含B网格点)赋予“预报”(水平平流和位移)插值。③、④、⑤ 步骤同于A-B流程。

4.3 Coons双三次曲面与Hermite双三次曲面片

 图 1 球面A网格上的3种拓扑矩形 (a)一般形, (b)临极点形, (c)极点形. P00、P01、P10、P11为变量值 Figure 1 Three topological rectangular figures of the spherical latitude-longitude grid. (a) general one, (b) one close to the pole, (c) one adjacent to the pole. In Fig. 1, P00, P01, P10 and P11 denote 4 variable values

4.4 求上游点

 $\frac{{p_s^{t + \Delta t} - p_s^t}}{{\Delta t}} = - \frac{{{{\tilde u}^\lambda } + {{\left( {\tilde v\cos \varphi } \right)}^\varphi }}}{{{r_0}\cos \varphi }} = - {{\tilde u}^x} - {{\tilde v}^y} + \frac{{\tilde v\tan \varphi }}{{{r_0}}},$ (62)

 $\int_A {\frac{{p_s^{t + \Delta t} - p_s^t}}{{\Delta t}}{\rm{d}}A} = - \int_A {\left[ {{{\tilde u}^\lambda } + {{\left( {\tilde v\cos \varphi } \right)}^\varphi }} \right]{r_0}{\rm{d}}\lambda {\rm{d}}\varphi } = 0,$ (63)

 $p_s^{t + \Delta t} = p_s^t + \Delta {{\tilde p}_s}\;\;\;\;;\Delta {{\tilde p}_s} = - \Delta t\left( {{{\tilde u}^x} + {{\tilde v}^y} - \frac{{\tilde v\tan \varphi }}{{{r_0}}}} \right),$ (64)

 $\left\{ \begin{array}{l} {\left( {\Delta {{\tilde p}_s}} \right)_N} = - \Delta t\left( {\tilde v_{N\left( {3{\rm{\pi /2}}} \right)}^y + \tilde v_{N\left( 0 \right)}^y} \right)\\ {\left( {\Delta {{\tilde p}_s}} \right)_S} = - \Delta t\left( {\tilde v_{S\left( {{\rm{\pi /2}}} \right)}^y + \tilde v_{S\left( 0 \right)}^y} \right) \end{array} \right..$ (65)
7 积分试验与分析

7.1 积分结果分析

A-B流程和A(B)-J流程积分结果非常接近(2008年1月10日积分结果见图 3~4, 2002年8月19日积分结果图略), 表明A-B流程积分过程(暂)未发生2倍格距波传播与增长, 即做A网格双三次曲面变换, 求取B网格上游点和J网格上游点在A网格上的水平坐标, 两者基本一致。

 图 3 2008年1月10日14:00初始全球地面(海平面)气压场(a), 以及以“lnp”为预报变量, 采用不同流程的“准守恒”与“守恒”积分10天的全球地面(海平面)气压场(单位: hPa) (a)初值场, (b) A(B)-J流程“静力质量守恒”格式, (c) A(B)-J流程“准守恒”格式, (d) A-B流程“静力质量守恒”格式 Figure 3 Initial pressure field at the global sea-level, conserved and quasi-conserved 10 days forecasted global sea-level pressure field applied "lnp" as forecast variable with different procedure at 14:00 on 10 January 2008. Unit: hPa. (a) initial pressure field, (b) conserved forecasted with A(B)-J procedure, (c) quasi-conserved forecasted with A(B)-J procedure, (d) conserved forecasted with A-B procedure
 图 4 2008年1月10日14:00初始全球500 hPa高度场(实线, 单位: dagpm)和温度场(虚线, 单位: K), 以及以“lnp”为预报变量, 采用不同流程的“准守恒”与“守恒”积分10天全球500 hPa高度场(实线, 单位: dagpm)和温度场(虚线, 单位: K) (a)初值场, (b) A(B)-J流程“静力质量守恒”格式, (c) A(B)-J流程“准守恒”格式, (d) A-B流程“静力质量守恒”格式 Figure 4 Initial global geopotential height field (solid line, unit: dagpm) and temperature field (dotted line, unit: K), conserved and quasi-conserved 10 days forecasted global geopotential height field (solid line, unit: dagpm) and temperature field (dotted line, unit: K) on 500 hPa at 14:00 on 10 January 2008. (a) initial pressure field, (b) conserved forecasted with A(B)-J procedure, (c) quasi-conserved forecasted with A(B)-J procedure, (d) conserved forecasted with A-B procedure

A(B)-J流程积分结果还表明, “准守恒”与“守恒”积分都考虑了全球三维静力平流, 两者仅地面气压预报方程有“些微”差别, 前者地面气压变率是对各层“质量散度”作气柱积分, 后者是对各层“质通量”作气柱积分之后求其平均散度, 但两者积分10天地面气压场和500 hPa高度场及温度场仍然非常接近(对比图 3图 4), 表明两者之间可能仅存在微弱的计算累积误差。但“准守恒”积分试验全球大气质量确有微小波动, 总趋势是地面平均气压从2008年1月10日初值为101133. 5 Pa(2002年8月19日初值为101140. 2 Pa)到积分10天的101129. 6 Pa(101139. 4 Pa), 且2008年1月10日地面平均气压最大波动振幅为101133. 5-101126. 7=6. 8 Pa(2002年8月19日为101140. 2-101135. 7=4. 5 Pa), 而“守恒”积分试验, 在绝对误差正负0. 1 Pa范围内, 全球大气质量无变化。另外, A-B流程“准守恒”积分试验地面平均气压有与A(B)-J流程相似的变化, 而其“守恒”积分试验地面平均气压也基本不变。

7.2 “失败”积分结果分析

 图 5 以“p”为预报变量, 采用A-B流程“静力质量守恒”格式积分6天的全球地面(海平面)气压场(单位: hPa) Figure 5 Conserved 6 days forecasted global sea-level pressure field applied "p" as forecast variable with A-B procedure. Unit: hPa
 图 6 以“p”为预报变量, 采用A-B流程“静力质量守恒”格式积分6天的500 hPa高度场(实线, 单位: dagpm)和温度场(虚线, 单位: K) Figure 6 Conserved 6 days forecasted global geopotential height field (solidline, unit: dagpm) and temperature field (dotted line, unit: K) at 500 hPa applied "p" as forecast variable with A-B procedure

8 结论

(1) 样条格式是关于全部格点的二阶可导非线性格式, 但样条格式线性部分是二阶中央差。本文在给定简谐波真解条件下, 推导证明二阶中央差比一阶中央差空间截断误差及相速和群速误差均减少一倍, 但两种格式相速和群速均与波数有关, 将真解非频散波拟合成为频散波, 且两种格式均为“群速慢于相速慢于真实气流”(能量逆向频散快于位相逆向传播), 波长越短, 误差越大, 这使得在扰动上游产生新的波动、或使得原有的波动增强, 出现上游/下游(向下游传播)效应, 这一方面对于预报天气系统移动、新生很重要, 另一方面也可能导致出现虚假“寄生(短)波”, 寄生波在基本气流上游逐渐形成, 并向下游传播, 寄生波发展可能造成某一变量场的计算紊乱, 但非线性样条格式的真实情况有待进一步(做理想场试验)研究。

(2) 与谱模式数学基础相当, 三次样条函数具备对于原函数及其1阶、2阶导数“收敛性”和对于原函数二阶导数最佳逼近“最优性”及“周期性”, 适合做全球/区域统一数值模式, 从而对应于谱模式动力框架核心, “高斯网格二维谱变换半隐式-半拉格朗日积分方案”, 本文提出样条模式动力框架核心, “准均匀经纬网格三次样条函数变换显式-准拉格朗日积分方案”。

(3) 给出准拉格朗日样条格式预报方程通式, 做“水平双三次曲面+垂直三次样条”拟合以求二阶时空精度上游点运动路径, 将广义牛顿力整体作为“匀加速”变率做风场预报, 以求一个时间步长水平位移以及“平均”散度场, 当作“变加速”变率做增压与增温(及水汽相变)预报, 从而实现全球静力守恒准均匀经纬网格样条格式数值模式, (而采用“三次曲体”拟合求上游点可做非静力数值模式), 这时, 气块的“三次运动路径”包含了“斜率”平流、“曲率”弯流和“挠率”扭流。风场预报(只)可以采用“匀加速”方案, 以避免采用半隐式-半拉格朗日积分方案, 进而避免求解Helmholtz方程, 而压、温(湿)场预报可采用计算精度较高的“变加速”方案, 因其变率就是“匀加速”方案已做预报的风场(位移)散度, 故整个“显式+隐式”方案, 本质上仍只是显式方案。

(4) 在球面经纬网格(A网格)基础上, 设计两种基本的准均匀经纬网格: B网格和J网格, 其中B网格点与A网格点重合, 但不可避免存在2倍格距过渡界线, 容易造成虚假的寄生(短)波传播与增长; J网格避免出现2倍格距过渡界线, 且J网格点一般不与A网格点重合。理论上还可设计多种准均匀经纬网格, 例如, 通过做“A网格-高斯网格”三次样条函数变换, 并通过做理想积分试验, 将来可与谱模式动力框架去比较计算精度。

(5) 做“A-B”网格三次样条函数变换积分试验, 若对对数气压(lnp)场做样条格式变换, 则没有发生虚假寄生波, 而若对气压(p)场做样条格式变换, 则出现明显的虚假寄生波传播与增长模拟个例, 发生这种现象, 是因为在A-B/A(B)-J/A-J网格上, lnp场与p场样条格式拟合斜率、曲率和挠率存在差别(可以肯定后者斜率较大), 其中A-B网格存在2倍格距过渡界线, 其一侧截断误差及相速、群速误差为另一侧的两倍, 从而促成了“失败”的模拟结果, 解决积分不稳定的方法有:在各种网格上做变量场平滑或样条格式光顺, 但平滑是“以精度换稳定”, 故将来应研究“准均匀(变)经纬网格”。

(6) 积分试验证明, 可设计出全球静力守恒样条格式准拉格朗日积分方案, 其做法是引入高度坐标静力连续方程, 其算法是直接对静力水平平流“位移”求样条格式散度场。

(7) 将来全球样条模式宜引入具有定常斜率、曲率和挠率的双三次曲面地形与地形追随高度坐标。

 Layton A T. 2002. Cubic spline collocation method for the shallow water equations on the sphere[J]. Journal of Computational Physics, 179(2): 578–592. DOI:10.1006/jcph.2002.7075 Bubnová R, Hello G, Bénard P, et al. 1995. Integration of the fully elastic equations cast in the hydrostatic pressure terrain-following coordinate in the framework of the ARPEGE/Aladin NWP system[J]. Mon Wea Rev, 123(2): 515–535. DOI:10.1175/1520-0493(1995)123<0515:IOTFEE>2.0.CO;2 Côté J, Gravel S, Méthot A, et al. 1998a. The operational CMC-MRB global environmental multiscale (GEM) model. Part Ⅰ:Design considerations and formulation[J]. Mon Wea Rev, 126(6): 1373–1395. DOI:10.1175/1520-0493(1998)126<1373:TOCMGE>2.0.CO;2 Côté J, Desmarais J G, Gravel S, et al. 1998b. The operational CMC|MRB global environmental multiscale (GEM) model. Part Ⅱ:results[J]. Mon Wea Rev, 126(6): 1397–1418. DOI:10.1175/1520-0493(1998)126<1397:TOCMGE>2.0.CO;2 Cullen M J P. 1990. A test of a semi-implicit integration technique for a fully compressible non-hydrostatic model[J]. Quart J Roy Meteor Soc, 116(495): 1253–1258. DOI:10.1002/(ISSN)1477-870X Cullen M J P, Davies T, Mawson M H, et al. 1997. An overview of numerical methods for the next generation UK NWP and climate model[J]. Atmos-Ocean, 35(supp1): 425–444. Ferguson J. 1964. Multivariable curve interpolation[J]. Journal of the ACM (JACM), 11(2): 221–228. DOI:10.1145/321217.321225 Gal-Chen T, Somerville R C J. 1975. On the use of a coordinate transformation for the solution of the Navier-Stokes equations[J]. Journal of Computational Physics, 17(2): 209–228. DOI:10.1016/0021-9991(75)90037-6 Li X, Chen D, Peng X, et al. 2006. Implementation of the semi-Lagrangian advection scheme on a quasi-uniform overset grid on a sphere[J]. Adv Atmos Sci, 23(5): 792–801. DOI:10.1007/s00376-006-0792-9 McDonald A, Bates J R. 1989. Semi-Lagrangian integration of a gridpoint shallow water model on the sphere[J]. Mon Wea Rev, 117(1): 130–137. DOI:10.1175/1520-0493(1989)117<0130:SLIOAG>2.0.CO;2 Ritchie H, Beaudoin C. 1994. Approximations and sensitivity experiments with a baroclinic semi-Lagrangian spectral model[J]. Mon Wea Rev, 122(10): 2391–2399. DOI:10.1175/1520-0493(1994)122<2391:AASEWA>2.0.CO;2 陈德辉, 胡志晋, 徐大海, 等. 2004. CAMS大气数值预报模式系统研究[M]. 北京: 气象出版社, 190. Chen Dehui, Hu Zhijing, Xu Dahai, et al. 2004. Research on the system of atmospheric numerical prediction model (CAMS)[M]. Beijing: China Meteorological Press, 190. 辜旭赞. 2003. 一模式大气参考状态的科学计算与特征分析[J]. 高原气象, 22(6): 608–612. Gu Xuzan. 2003. Compution and characteristic of a model atmospheric reference state[J]. Plateau Meteor, 22(6): 608–612. 辜旭赞. 2010. "准拉格朗日法"和"欧拉法"数学一致性与三次插值函数算法时间积分方案[J]. 高原气象, 29(3): 655–661. Gu Xuzan. 2010. The mathematical consistency of quasi-Lagrangian and Eulerian time integratifmn schemes with fitting ciibic interpolation functifmn[J]. Plateau Meteor, 29(3): 655–661. 辜旭赞. 2011. 一种"准拉格朗日法"和"欧拉法"统一算法时间积分方案[J]. 气象学报, 69(3): 440–446. DOI:10.11676/qxxb2011.038 Gu Xuzan. 2011. A new quasi-Lagrangian time integration scheme with the interpolation of fitting bicubic surface[J]. Acta Meteor Sinica, 69(3): 440–446. DOI:10.11676/qxxb2011.038 李江浩, 彭新东. 2013. 阴阳网格上质量守恒计算性能分析[J]. 大气科学, 37(4): 852–862. DOI:10.3878/j.issn.1006-9895.2012.12060 Li Jianghao, Peng Xindong. 2013. Analysis of computational performance of conservative constraint on the Yin-Yang Grid[J]. Chinese J Atmos Sci, 37(4): 852–862. DOI:10.3878/j.issn.1006-9895.2012.12060 廖洞贤. 1999. 大气数值模式的设计[M]. 北京: 气象出版社, 291. Liao Dongxian. 1999. Design of atmospheric numerical model[M]. Beijing: China Meteorological Press, 291. 吕美仲, 欧阳子济, 翟子航. 1983. 正压原始方程半隐式样条格式[J]. 气象科学, 4(2): 1–11. Lü Meizhong, Ouyang Ziji, Zhai Zihang. 1983. The semi-implicit spline scheme of the barotropic primitive equation[J]. J Meteor Sci, 4(2): 1–11. 奚梅成. 1995. 数值分析方法[M]. 合肥: 中国科学技术大学出版社, 377. Xi Meicheng. 1995. Numerical analysis method[M]. Hefei: University of Science & Technology China Press, 377. 薛纪善, 陈德辉. 2008. 数值预报系统GRAPES的科学设计与应用[M]. 北京: 科学出版社, 383. Xue Jishan, Chen Dehui. 2008. Scientific design and application of the numerical prediction system (GRAPES)[M]. Beijing: Science Press, 383. 张玉玲, 吴辉碇, 王晓林. 1986. 数值天气预报[M]. 北京: 科学出版社, 472. Zhang Yulin, Wu Huizhan, Wang Xiaolin. 1986. Numerical weather Prediction[M]. Beijing: Science Press, 472. 周毅, 刘宇迪, 桂祁军, 等. 2002. 现代数值天气预报[M]. 北京: 气象出版社, 220. Zhou Yi, Liu Yudi, Gui Qijun, et al. 2002. Modern numerical weather prediction[M]. Beijing: China Meteorological Press, 220.
A Quasi-Lagrangian Integration of Conservation of Atmospheric Mass with Unify Scheme of Cubic Spline Function Transformating on Quasi-uniform Latitude-longitude Grid and Its Integration Cases
GU Xuzan1 , ZHAO Jun2 , TANG Yonglan1
1. Hubei Key Laboratory for Heavy Rain Monitoring and Warning Research, Institute of Heavy Rain, China Meteorological Administration, Wuhan 430074, China;
2. Institute of Ocean Science and Engineering, National University of Defence Technology, Changsha 410073, China
Abstract: The spline format is a no-linear, second-order derivative one, its linear segment is that of the second-order central difference. In this paper, we give a derivation proof of that the space truncation error, phase velocity and group velocity errors of the second-order center differential is halved that of the first-order center differential under a hypothesis of genuine solution of simple harmonic wave. So, we draw lessons from the idea of the dynamic core of spectral model, the semi-implic semi-Lagrangian integration scheme with 2D spectral spherical harmonic function transform on the Gussion grid, to introduce a new explicit quasi-Lagrangian integration scheme with cubic spline function transform on guasi-uniform latitude-longitude grid (called "spline model"). Adopting original atmospheric equations of motion, which includes that in the North Pole and the Sorth Pole, a general forecast equation of spline format of space-time second-order differential remainder is derived, then obtain the hydrostatic pressure and temperature forecast eqations of conservation of the atmospheric mass. Based on uniform latitude-longitude grid, we harmonize two quasi-uniform ones, which must be quasi-uniform latitude space, and on which cubic spline function transformation (transformation=fitting+interpolation) must be done for variables of pressure, temperature, moisture, winds and general Newtonian force acting to unit air mass on rotating earth (acceleration), which made all of them second-order derivative, to solve the track of an upstream point, but the upstream air parcel goes alone just "cubic path" of fitting their slopes, curvatures and torsions of the variable fields to "bicubic surface in horizontal + cubic spline in vertical". It is with a path of uniform acceleration motion to forecast wind field, and with fitting splines to the paths of the 3D hydrostatic advection and getting its implicit average divergence in one time step, to forecast increments of pressure and temperature fields in the adiabatic process. We give two integration cases that testify to the dynamic core of global spline model.
Key Words: Cubic spline function (spline format)    Quasi-uniform latitude-longitude grid    Space-time second-order discrete forecast equation    Quasi-Lagrangian integration scheme    Iintegration cases