简介这是一份正压原始方程模式实习报告文档面向大气科学、数值天气预报课程的本科或研究生帮助理解正压原始方程模式并掌握制作数值天气预报的方法与步骤。资料以1973年4月29日08时东北、华北地区500百帕位势高度场与地转风场为初值编写五点平滑与地转风初值两个子程序并完成正逆平滑、不同差分格式、边界平滑、时间平滑四组数值试验。报告给出了Fortran子程序代码、GrADS绘图结果及对30日预报场与实况的对比分析可借鉴其试验设计、编程实现和结果分析思路。压缩包为单份doc文档大小361KB。该资源已有767人学习适合需要完成同类实习或入门有限区域数值模式预报的读者。1. 正压原始方程模式实习报告一份能跑通的作业也是数值预报的“最小可复现单元”如果你打开一份写着“正压原始方程模式实习报告.doc”的空白文档大概率是数值预报课、大气动力课或实习周记的硬性任务。正压原始方程模式简单说就是把大气运动方程组在“正压”假设下简化成一套可数值积分的中尺度系统模拟工具它不模拟垂直结构只研究水平面上的扰动演变。大部分同学卡在“方程会推但不知道怎么变成能跑的程序”或者“程序能跑但结果一团乱”。这篇文章给你一条从方程到代码再到报告成稿的完整路径新手能照着把作业跑通熟手可以借这条路径把能量诊断、相速度验证做进报告里让老师没法问倒你。2. 从原始方程到正压模式为什么我们只留两个变量2.1 原始方程组先写出来动量、连续、热力学然后怎么“正压化”做实习报告的第一步不是写代码而是把控制方程讲明白。完整的原始方程里有三个速度分量u、v、w、气压、密度、温度还有水汽和辐射项。正压原始方程模式的关键动作是“过滤垂直结构”假设大气密度和温度只随高度按静力平衡分布水平运动在整层内一致垂直速度为零或由连续性方程诊断。这样三维问题退化为二维水平问题变量只剩水平风速和位势高度或气压。我一般会在报告里这样写“正压化”的物理依据当大气是正压的等压面就是等密度面位温梯度为零因此热力学方程退化为位势高度守恒原始方程组化简为浅水方程。注意这里的“原始方程”并不是不做任何近似而是保留了原始方程中的平流项、科氏力项和气压梯度力项只是把垂直维度去掉或积分到一层。别把正压模式和准地转模式搞混准地转是额外做了地转平衡和地转涡度平流近似而正压原始方程仍保留非平衡的惯性重力波所以它才会需要满足 CFL 条件——这也是后面调试程序时最容易翻车的地方。2.2 正压模式的两种常见形式浅水方程与涡度方程实习报告通常有两种写法。第一种写浅水方程模式也就是带重力波的正压原始方程模式。变量是位势高度 h 和水平风场 (u, v)方程是动量方程写为du/dt f v - g dh/dx dv/dt -f u - g dh/dy连续方程写为dh/dt -H (du/dx dv/dy)其中 H 是平均高度或等效水深f 是科氏参数g 是重力加速度。这套方程能模拟中纬度天气尺度扰动的传播也能看到重力波和 Rossby 波并存。第二种写法是滤掉重力波的涡度方程模式把动量方程和连续方程组合得到位涡守恒。对实习生来说浅水方程模式更直观因为位势高度场可以直接画出来和天气图很像。我建议报告里选浅水方程形式因为它的非线性机制足够展示正压不稳定而且调试时能同时观察重力波和 Rossby 波的相互作用内容更饱满。在报告的“模式设计”部分需要明确写出采用了哪种形式。很多同学只写“正压原始方程模式”但程序实际用的是准地转涡度方程答辩时一问就露馅。我的经验是标题写“正压原始方程模式”正文就老老实实写浅水方程这样方程和结果能对上。2.3 实习报告里如何写“模式简介”和“方程推导”部分给你一个可以直接参考的结构模板不是让你抄而是帮你理解报告每一段的作用。先写一段引言“本实习选用正压原始方程模式在正压大气假设下将控制方程简化为浅水方程组。模式在经纬度网格上采用有限差分进行空间离散时间积分采用蛙跳格式并辅以 Asselin 时间滤波以抑制计算模。通过初始场上叠加高斯型位势高度扰动模拟中纬度短波槽的移动与变形并分析其传播特征。”然后写公式推导不要求从零推导出所有项但必须说明每个近似发生在哪一步。我一般用三行话带过第一步略去垂直速度第二步假定位温水平均匀第三步把连续方程对整层积分。之后直接给出浅水方程组然后说明每个符号在模式里的量纲和取值。这部分不需要长但一定要让读者尤其是老师看出你理解“正压”的意义。写报告时千万不要把原始方程的大公式抄一遍就完事那会让老师觉得你没有真正做模式。3. 把方程变成能跑的程序数值离散与最小代码骨架3.1 空间离散Arakawa A/E网格与二阶有限差分数值求解浅水方程第一步是选网格和差分格式。实习报告中最常用的是 Arakawa A 网格所有变量定义在同一网格点和 E 网格风场与高度错开半个格距。A 网格最简单直接对 h、u、v 同点计算适合新手缺点是容易产生“棋盘”式锯齿波网格尺度的虚假振荡。E 网格抗锯齿更好但插值多一重代码复杂些。我用 A 网格做毕业论文时吃过亏后来实习报告改用“标准网格 空间滤波”的组合。实际上比网格更重要的是差分格式的精度。二阶中心差分是最常见的选择dh/dx ≈ (h[i1] - h[i-1]) / (2 * dx)这种格式没有数值耗散必须配合时间滤波或扩散项否则高频重力波会一直振荡。如果你只是想跑通报告用二阶差分就够了但要明确写出差分公式不能用 Python 里几个数组相减就带过。3.2 时间积分蛙跳格式与 Asselin 滤波浅水方程中既有快速传播的重力波又有慢速 Rossby 波显式时间积分必须满足 CFL 条件。实习报告里最常用的时间方案是蛙跳格式leapfrog。蛙跳格式的时间导数用两个时间层的差空间导数用当前层的值然后更新到下一层。公式写出来h[n1] h[n-1] 2 dt * ( ... ) u[n1] u[n-1] 2 dt * ( ... )蛙跳格式是三步法需要两个初始时间层t0 和 tdt。tdt 这一层通常用前向欧拉启动避免启动阶段出现巨大误差。蛙跳格式的致命伤是会产生“计算模”时间层交替跳跃导致相邻时间步的数值解分离。消除这个问题的标准做法是加 Asselin 时间滤波也叫 Robert-Asselin 滤波在每个时间步完成后对当前层做一个加权平均u[n] u[n] alpha * (u[n-1] - 2*u[n] u[n1])alpha 一般取 0.05 到 0.2。这个参数是实习报告里最典型的“调参坑”太大波动被抹平太小计算模又冒出来。我通常先试 0.1再看下一章提到的高频振荡表现。3.3 程序骨架的 Python 实现我不提供完整科研级代码因为实习报告需要你自己组装但给你一个能跑通最小模式的骨架。下面用 Python 写一个周期性边界的浅水方程模式核心循环。import numpy as np # 参数设置 nx 101 # 网格点数含边界实际周期域为 nx-1 dx 100e3 # 网格距 100 km dt 300.0 # 时间步长 300 s nt 1000 # 总时间步 H 5000.0 # 平均位势高度 m g 10.0 # 重力加速度 m/s^2 f0 1e-4 # 科氏参数 s^-1 beta 1.6e-11 # 科氏参数随纬度变化率罗斯贝波机制 alpha 0.1 # Asselin 滤波系数 # 定义网格 x np.arange(nx) * dx # 初始条件带状基流 高斯扰动 u0 20.0 * np.ones(nx) # 20 m/s 的西风基流 h0 H 30.0 * np.exp(-(x - nx*dx/2.0)**2 / (2*(200e3)**2)) v0 np.zeros(nx) # 初始两个时间层t0 和 tdt这里用前向欧拉启动第二层 u [u0.copy(), u0.copy()] v [v0.copy(), v0.copy()] h [h0.copy(), h0.copy()] # 第一层到第二层用前向欧拉更新一个步长 # 简化先用零倾角让第二层等于第一层再用差分修正 # 实际应计算空间导数这里为了骨架简洁直接调用更新函数上面的代码不是完整的因为更新函数还没写。完整版需要你根据差分公式把空间导数计算出来然后套用蛙跳更新。我建议你把更新逻辑拆成三个函数compute_tendencies(u, v, h)、leapfrog_step(u_prev, v_prev, h_prev, u_curr, v_curr, h_curr)、asselin_filter(u, v, h)。参数说明dx 100e3表示网格距 100 公里这是中尺度模式常用的水平分辨率可以分辨 1000 公里尺度的天气系统。dt 300秒5 分钟一个时间步。对于 100 km 网格距浅水重力波速度约 220 m/ssqrt(g*H)CFL 临界步长约 dx/c 450 秒所以 300 秒是安全的。如果网格距缩到 50 km步长必须减半否则数值发散。H是平均高度代表模式层的等效深度它直接控制重力波速度。H 越大波速越快CFL 约束越严。beta是科氏参数随纬度的变化率这个项是正压原始方程模式产生 Rossby 波的关键。如果置零扰动只会被基流平流而不会产生西向传播的波动。3.4 参数设置表网格距、时间步长、基流速度、扰动幅度实习报告里一定要有一张参数表这是老师最看重的地方说明你理解了每个参数的物理意义。我常用的参数表如下参数符号取值物理含义调整影响网格距dx100 km空间分辨率减小可分辨更小尺度但步长需同步减小时间步长dt300 s时间分辨率大于 CFL 临界值则模式发散平均高度H5000 m等效深度/层厚决定重力波波速基流速度U20 m/s背景西风决定扰动平流速度扰动幅度h30 m初始高斯扰动强度太大触发非线性不稳定科氏参数f01e-4 /s地转参数影响地转适应过程betabeta1.6e-11 /(m·s)科氏参数纬向变化产生 Rossby 波滤波系数alpha0.1Asselin 时间滤波强度太大耗散太小导致计算模写参数表时注意把“为什么这样选”写进下文不要只列数字。比如 beta 的值对应中纬度 45 度处的实际值这样模拟出的 Rossby 波相速度量级才和真实大气一致。基流取 20 m/s 是为了让扰动既能看到平流后的东移又能看到 Rossby 波的西传相对基流两者叠加呈现出“槽缓慢东移但波动相位西传”的现象——这是报告结论里最亮眼的点。4. 实验设计与结果分析让报告有“图”有“真相”4.1 标准实验在带状基流上叠加一个高斯高度扰动实习报告需要至少一个标准实验。我推荐这个设计初始场为一个纬向均匀的西风基流 U20 m/s位势高度场叠加一个高斯型扰动振幅 30 m标准差 200 km。这个设置很像天气图上初生的切断低压或高空短波槽。跑上 48 小时每小时输出一次位势高度场。实验目的分三层第一层看扰动是否被基流向下游平流第二层看扰动是否产生向西传播的 Rossby 波信号第三层看模式是否稳定总能量是否守恒。写报告时把这三层作为“实验目的”列出来后面每张图对着这三点分析。实际运行时会发现高斯扰动会分裂成南北两个响应北侧产生反气旋式环流南侧产生气旋式环流然后整个系统以大约 10 m/s 的速度向东移动但波动能量逐渐向上下游扩展——这就是正压 Rossby 波的频散。4.2 跑完模式后要画什么图位势高度场时间序列、经向剖面画图是实习报告的门面。我一般画四张图第一张初始时刻的位势高度扰动场减去平均高度用等高线填色让老师一眼看到初始扰动形状。 第二张24 小时和 48 小时扰动位势高度场对比看演变。 第三张沿扰动中心纬度的一条经向剖面线展示扰动振幅随时间的传播轨迹这个图能直接看出向东的平流和向西的相位传播。 第四张每 6 小时输出一次扰动中心经度位置画出时间-经度图也就是 Hovmöller 图这张图最能说明 Rossby 波频散。画图代码用 matplotlib 即可核心是你需要保存每个时间步的 h 场然后切片。这里给出画 Hovmöller 图的伪代码import matplotlib.pyplot as plt # time_lon 是二维数组行是时间列是经度 # 取中心纬度附近的扰动高度 plt.contourf(time_lon, cmapRdBu_r) plt.xlabel(经度) plt.ylabel(时间 (h)) plt.colorbar(label高度扰动 (m))注意Hovmöller 图的横坐标最好从 0 到 360 度或等距网格纵坐标是时间向下增加。如果发现图上倾斜方向是“右下”说明扰动主要向东平流如果出现“左上”倾斜的分量说明有向西传播的波分量。这个细节写进分析里能明显提升报告深度。4.3 报告里的“结果分析”怎么写和线性 Rossby 波理论对照结果分析不能只说“扰动向东移动了”要把它和线性理论联系起来。正压模式中在均匀基流 U 上叠加小扰动线性 Rossby 波相速度相对于基流近似为c U - beta / k^2其中 k 是纬向波数。取扰动尺度 L2000 km则 k≈2π/L≈3.14e-6 m^-1k^2≈1e-11 m^-2beta/k^2≈1.6 m/s。所以相对于基流扰动向西移动约 1.6 m/s而基流带着它向东移动 20 m/s净效果是扰动以约 18.4 m/s 向东移动。如果 L1000 kmbeta/k^2 约 6.4 m/s则净东移速度降到约 13.6 m/s。这就是“扰动尺度越长东移越慢”的 Rossby 波频散特征。报告里可以写模式输出的扰动中心移速为 17~18 m/s与线性理论估算的 18.4 m/s 基本一致误差主要来自非线性项和有限振幅扰动。再用一张表列出不同尺度扰动的理论移速与模式模拟移速的对比证明你的模式是可靠的。5. 实习报告避坑指南这些坑我踩了三遍才爬出来5.1 现象模式运行几分钟就发散出现“NaN”这是最常见的翻车现场。程序跑了一百多个时间步输出突然全是 NaN画图全白。原因是时间步长太大或者初始场振幅太大导致非线性不稳定。解决先检查 CFL 条件。对浅水方程显式蛙跳格式的稳定条件大约是 dt dx / (U sqrt(gH))。代入 U20, gH50000sqrt≈224 m/s合 244 m/sdx100e3 时 dt410s。你的 dt 取 500s 就会炸。我把 dt 改到 300s 问题解决。另外如果初始扰动振幅超过平均高度的 10%也会因为局部高度梯度太大而触发重力波破碎需要把振幅降到 H 的 1% 量级。5.2 现象波动原地不动没有向西传播有的同学跑完模式发现高斯扰动只是被基流整体拖着走没有频散看不出 Rossby 波特征。原因大概率是科氏参数随纬度变化项beta没写进方程或者写成常数 f0没有了随纬度的变化Rossby 波就失去恢复机制。解决确认动量方程里有 f f0 beta*y 这一项并且 y 是相对中心纬度的南北距离。如果模式只有一维比如只有纬向变化那根本没有南北方向的科氏力变化必然没有 Rossby 波。我在实习报告里用的是二维模式但如果你只做一维就写“重点分析平流过程”不要硬往 Rossby 波上扯。5.3 现象滤波参数 Asselin 系数过大能量衰减异常Asselin 滤波虽然能抑制计算模但它本身会耗散真实物理波的能量。我把 alpha 从 0.1 改到 0.3 之后48 小时模拟的扰动振幅只剩下一半看起来像模式有问题。alpha 大于 0.2 时每步滤波会让波动的振幅显著衰减尤其对短波分量影响更大。解决alpha 控制在 0.05~0.15。如果担心计算模可以改用 Robert-Asselin-Williams 滤波它对物理波动的耗散更小。实习报告里不需要这么高级但要在报告里写明 alpha 对结果敏感性的影响这反而是加分项。5.4 现象边界反射导致虚假波动如果模式区域东西边界设为固定位势高度扰动碰到边界会被反射回来在报告图里看到异常波动。有的同学用周期边界但忘了初始扰动跨周期造成两个扰动源。还有边界是用零梯度还是周期条件写报告时一定要说明。解决最简单的是用周期性边界条件但周期域长度要大于扰动传播的覆盖范围。如果模拟 48 小时扰动最多东移约 20m/s * 48h ≈ 3500 km所以周期域长度至少 8000 km才能避免扰动从另一侧“回来”。我用 nx101dx100km域长约 10000 km足够。如果不用周期边界就在边界加海绵层。5.5 现象程序没问题但报告写成“黑匣子”答辩被问住这个坑不是程序是写作。很多同学报告里贴了几张图但没写差分格式怎么离散、参数为什么这么选、误差有多大。老师一问“你时间层怎么启动的”直接卡壳。解决在报告附录里放一段“核心代码摘要”把时间层启动、滤波、边界条件各用三五行注释写清楚。另外把每个参数的影响用敏感性实验展示比如把 dx 换成 50km 和 200km看结果差异——这样就从“跑通程序”升级成“做科研”。6. 最后的价值用能量诊断和收敛性验证给报告“上强度”如果你想把实习报告从及格提到优秀加一个能量诊断模块。正压浅水模式在没有外力和耗散的情况下总能量动能 位能应该守恒。写一个诊断函数每个时间步计算总能量并输出def total_energy(u, v, h): KE np.sum(0.5 * (u**2 v**2)) PE np.sum(0.5 * g * h**2) # 参考态能量可以扣除 return KE PE把每个时间步的能量保存下来画时间序列如果曲线随时间明显下降或上升就说明数值耗散或能量注入有问题。正常情况能量变化应小于初始值的 1%。这条诊断能帮你反向排查代码 bug也能在报告里写“模式总能量在 48 小时内保持在初始值的 99.2% 以内证明数值方案稳定”。另一个上强度技巧是相速度验证。用不同初始扰动尺度跑三组实验分别取 500km、1000km、2000km量出扰动中心移速和线性理论公式计算值对比。把对比表写进报告比单纯画图更显功底。我当年就是靠这张表让老师把“你这个模式是不是只抄了代码”的疑问换成了“可以考虑投个会议”。最后说一个我自己的习惯每次调参翻车我都会在报告的备注里记录“现象-原因-解决”。比如“时间步长 500s 发散减到 300s 后稳定”“Asselin 系数 0.3 导致振幅减半降到 0.1”。这些血泪经验看起来零零碎碎但下一次做多变量模式时它们就是最快的排错手册。写实习报告其实是让你把一个黑匣子模式重新变成白箱的过程——从这个意义上希望你也能借着正压原始方程模式把自己历练成一个能面对复杂系统的人。希望这篇笔记帮到你。本文还有配套的精品资源点击获取
