企业官网建设流程全解析

简介：这份正压原始方程模式实习报告面向大气科学、气象学及数值天气预报方向的本科生与研究生，用于课程实习、上机实验与报告撰写参考。报告以1973年4月29日08时东北、华北地区500百帕位势高度场及地转风场为初值，采用固定水平侧边界条件与二次守恒平流格式，制作未来24小时有限区域高度场和风场预报，并围绕五点平滑子程序、地转风初值子程序展开四个数值试验，涵盖正逆平滑对比、不同差分格式、边界平滑与时间平滑等内容，附有程序代码、预报图形及结果分析。资源包共1个doc文件，约361KB，结构完整、便于直接查阅与二次编辑。目前已有767人学习下载，适合需要理解正压原始方程模式流程、掌握数值天气预报一般方法与步骤、并训练Fortran编程与GrADS绘图的读者参考使用。

1. 正压原始方程模式实习报告：一份能跑通 24 小时预报的 Fortran 源码包

如果你手头正压着一份数值天气预报的课程实习，要求用 1973 年 4 月 29 日 08 时东北、华北 500 百帕位势高度场做初值，跑出未来 24 小时有限区域的高度场和风场，那你大概率绕不开这份正压原始方程模式实习报告。它不是纯理论推导，而是一套能直接编译运行的 Fortran 源码加实验记录：两个核心子程序（五点平滑、地转风初值）、四个数值试验（正平滑对比、差分格式对比、边界平滑、时间平滑），外加 GrADS 出图。适合正在做数值预报实习的学生、需要快速复现正压原始方程模式的从业者，以及想拿一份完整代码框架改造成自己实验的人。我拆过不少类似的实习包，这份的代码结构算清晰的，但坑也埋得实在。

2. 正压原始方程模式的核心：二次守恒平流格式与两个子程序

2.1 为什么选正压原始方程而不是斜压模式

正压原始方程模式建立在 500 百帕等压面上，假设大气是正压的、准地转的、准静力平衡的。这意味着密度只依赖气压，温度场和风场没有垂直切变，大气环流形势相对简单稳定。对于 24 小时有限区域预报来说，这个假设带来的误差在可接受范围内，尤其是影响中国的主要系统——贝加尔湖低压槽——预报效果反而不错。斜压模式当然更真实，但计算量大、初值要求高，实习阶段用正压模式能把数值计算、编程、上机操作的完整链路走一遍，性价比最高。

二次守恒平流格式是这套模式的核心。普通差分格式在长时间积分中会虚假制造或耗散能量，导致预报场失真。二次守恒格式保证动能和位涡拟能的守恒性，让 24 小时积分后场的基本形态不崩。代码里没有直接写“二次守恒”四个字，但平流项的差分构造和五点平滑的配合，实际上就是在维持这个守恒性。

2.2 五点平滑子程序：正平滑与正逆平滑的实现细节

五点平滑子程序ssip是这份代码里第一个要啃的硬骨头。它的作用是对格点场做空间滤波，抑制短波噪声。子程序接收数组a(m,n)、工作数组w(m,n)、平滑系数s、格点数m,n和开关k,l。核心逻辑是：对每个内点，用周围四个点和中心点的加权平均更新值。

subroutine ssip(a,w,s,m,n,k,l) dimension a(m,n),w(m,n) if(l==1) then do i=2,m-1 do j=2,n-1 w(i,j)=a(i,j)+s*(a(i-1,j)+a(i+1,j)+a(i,j-1)+a(i,j+1)-4*a(i,j))/4.0 end do end do do i=2,m-1 do j=2,n-1 a(i,j)=w(i,j) end do end do return else do i=2,m-1 do j=2,n-1 w(i,j)=a(i,j)+s*(a(i-1,j)+a(i+1,j)+a(i,j-1)+a(i,j+1)-4*a(i,j))/4.0 end do end do do i=2,m-1 do j=2,n-1 a(i,j)=w(i,j) end do end do do i=2,m-1 do j=2,n-1 w(i,j)=a(i,j)+(-s)*(a(i-1,j)+a(i+1,j)+a(i,j-1)+a(i,j+1)-4*a(i,j))/4.0 end do end do do i=2,m-1 do j=2,n-1 a(i,j)=w(i,j) end do end do endif return end

逻辑说明：l==1时只做一次正平滑，l取其他值时先正平滑再逆平滑。正平滑的公式是a + s*(四邻点之和 - 4a)/4，本质是拉普拉斯滤波。逆平滑把系数取反再算一遍，目的是在滤掉短波的同时尽量保留原始场的振幅。参数s一般取 0.5 左右，太大场会过度平滑，太小滤波效果不够。m,n是格点维度，循环从 2 到 m-1、2 到 n-1，边界点不参与平滑，这是固定侧边界条件的要求。

2.3 地转风初值子程序：差分格式与边界处理

地转风初值子程序cgw负责从位势高度场za反算地转风ua,va。公式用的是地转风关系：u = -(1/f) * ∂Φ/∂y，v = (1/f) * ∂Φ/∂x，代码里乘了rm（地图放大系数）和9.8（重力加速度），除以f（科氏参数）和d（格距）。

subroutine cgw(ua,va,za,rm,f,d,m,n) implicit none integer m,n,p,q,i,j real ua(m,n),va(m,n),za(m,n),rm(m,n),f(m,n),d p=m-1 q=n-1 do i=1,m ua(i,1)=-rm(i,1)*9.8*(za(i,2)-za(i,1))/(f(i,1)*d) ua(i,n)=-rm(i,n)*9.8*(za(i,n)-za(i,n-1))/(f(i,n)*d) do j=2,q ua(i,j)=-rm(i,j)*9.8*(za(i,j+1)-za(i,j-1))/(2.0*f(i,j)*d) end do end do do j=1,n va(1,j)=rm(1,j)*9.8*(za(2,j)-za(1,j))/(f(1,j)*d) va(m,j)=rm(m,j)*9.8*(za(m,j)-za(m-1,j))/(f(m,j)*d) do i=2,p va(i,j)=rm(i,j)*9.8*(za(i+1,j)-za(i-1,j))/(2.0*f(i,j)*d) end do end do end subroutine cgw

逻辑说明：ua是东西风分量，va是南北风分量。边界点用单侧差分（前差或后差），内点用中央差分。注意ua(i,1)和ua(i,n)用的是za(i,2)-za(i,1)和za(i,n)-za(i,n-1)，这是沿 y 方向的一阶差分；va(1,j)和va(m,j)用的是沿 x 方向的一阶差分。内点ua(i,j)用za(i,j+1)-za(i,j-1)除以2d，va(i,j)用za(i+1,j)-za(i-1,j)除以2d。参数d是格距，f随纬度变化，rm是地图投影放大系数。差分格式的选择直接影响初值场的平滑程度，前差和后差在边界附近会产生不同的风场结构，这就是后面数值试验要对比的内容。

3. 四个数值试验怎么跑：从编译到出图的完整流程

3.1 编译环境与依赖准备

这份代码是 Fortran 77 风格，用gfortran编译最省事。GrADS 用来画图，需要提前装好并配置GADDIR环境变量。我一般会建一个工作目录，把源码、初值数据、GrADS 脚本分开放。

mkdir -p ~/barotropic_exp/{src,data,grads} cd ~/barotropic_exp/src # 把主程序和两个子程序保存为 model.f gfortran -o model model.f -O2 -fcheck=all

-O2开优化，-fcheck=all打开运行时检查，数组越界和未初始化变量会直接报错，比事后查玄学问题省事得多。编译通过后，把初值数据按代码要求的格式写成文本文件，一般是m行n列的位势高度场。

3.2 试验一：正平滑与正逆平滑对比

这个试验改的是ssip子程序的l参数。l=1只做正平滑，l=2做正逆平滑。跑两次，分别输出 24 小时预报场，用 GrADS 叠在一起看差异。

# 正平滑 ./model < input_l1.txt > output_l1.txt # 正逆平滑 ./model < input_l2.txt > output_l2.txt

正平滑会持续耗散短波，场变得过于光滑，低压中心强度可能被削弱。正逆平滑在滤掉噪声的同时把振幅拉回来一部分，预报场更接近实况。但逆平滑不是万能的，如果s太大，正逆平滑也会引入新的计算噪声。我一般先用s=0.5跑一组，看预报场和实况的偏差再微调。

3.3 试验二：地转风差分格式对比

这个试验改的是cgw子程序里边界点和内点的差分方式。代码里已经写了前差和后差两种边界处理，内点固定中央差分。你可以把边界点的单侧差分改成中央差分（需要虚拟点），或者把内点的中央差分改成前差/后差，观察初值风场的变化。

! 把内点中央差分改成前差 ua(i,j)=-rm(i,j)*9.8*(za(i,j+1)-za(i,j))/(f(i,j)*d) va(i,j)=rm(i,j)*9.8*(za(i+1,j)-za(i,j))/(f(i,j)*d)

前差会让风场整体偏移，后差偏移方向相反，中央差分最平衡。边界点用单侧差分是必须的，因为边界外没有数据。这个试验的价值在于让你直观看到差分格式对初值场的敏感性——初值差一点，24 小时后预报场可能差很多。

3.4 试验三与试验四：边界平滑和时间平滑的开关

边界平滑控制的是ssip子程序里循环的范围。默认只平滑内点，边界点不动。如果把循环改成do i=1,m和do j=1,n，边界点也参与平滑，但固定侧边界条件要求边界值不变，所以平滑后要把边界值还原。时间平滑是在时间积分过程中每隔几步做一次平滑，抑制高频振荡。代码里没有直接给时间平滑的子程序，常见做法是在时间步循环里插入ssip调用，每 6 步或 12 步平滑一次。

! 时间步循环内插入平滑 if (mod(nt, 6) == 0) then call ssip(h, w, 0.5, m, n, 1, 1) end if

nt是时间步计数，每 6 步做一次正平滑。不做时间平滑时，24 小时积分后可能出现明显的短波噪声，高度场等值线变得锯齿状。做了时间平滑，场更干净，但过度平滑会抹掉真实的中尺度系统。

4. 避坑与排查：这份代码跑不起来时先看这几条

4.1 数组越界导致段错误

现象：编译通过，运行几秒后报Segmentation fault，没有任何其他提示。原因：ssip和cgw子程序里循环边界写的是2,m-1和2,n-1，但如果主程序传入的m,n和数组实际维度不一致，或者初值数据行列数读错，就会越界。解决：在编译时加-fcheck=all，运行时会打印越界的具体位置。同时检查读初值的read语句，确保m,n和数组声明一致。

4.2 平滑系数 s 取太大导致场崩

现象：预报场在几小时后出现明显的棋盘格噪声，等值线密集到无法辨认。原因：s超过 0.5 后，五点平滑的滤波特性变差，正逆平滑的逆平滑项会放大高频噪声。解决：s从 0.3 开始试，逐步加到 0.5，观察预报场形态。如果噪声已经出现，把s降到 0.2 再跑一次，对比噪声是否消失。

4.3 地转风初值在边界附近异常大

现象：初值风场在区域边界出现风速极值，远超合理范围。原因：边界点用单侧差分，如果边界附近的位势高度梯度本身很大，除以d后风速会被放大。另外f在低纬度接近零，如果区域包含低纬度，1/f会爆炸。解决：检查区域范围，确保f不在零附近。边界点的差分可以改用二阶精度单侧差分，或者对边界风场做一次额外平滑。

4.4 GrADS 出图时变量不匹配

现象：GrADS 能打开数据文件，但d h显示全零或报错。原因：Fortran 写二进制数据时的记录长度和 GrADS 控制文件.ctl里的undef、xdef、ydef不匹配。解决：确认 Fortran 输出是直接访问还是顺序访问，GrADS 读直接访问二进制需要.ctl里写options big_endian或little_endian。先用ncdump或od看数据文件的字节数，反推记录长度。

4.5 时间积分步长过大导致计算不稳定

现象：积分到第 10 步左右，高度场数值突然增大几个量级。原因：正压原始方程的 CFL 条件要求时间步长小于格距除以最大风速。如果d是 300 公里，风速 50 米/秒，时间步长不能超过 6000 秒。解决：把时间步长降到 300 秒或 600 秒，重新积分。如果必须用大步长，需要改用半隐式格式，但这份代码是显式格式，只能缩步长。

5. 进阶技巧：用 GrADS 脚本自动对比四组试验

跑完四个试验后，手动一张张看图效率太低。我一般写一个 GrADS 脚本，把四组预报场画在一张图上，直接看差异。

'reinit' 'sdfopen output_l1.ctl' 'sdfopen output_l2.ctl' 'sdfopen output_l3.ctl' 'sdfopen output_l4.ctl' 'set lon 100 140' 'set lat 30 55' 'set lev 500' 'set gxout contour' 'set ccolor 1' 'd h.1' 'set ccolor 2' 'd h.2' 'set ccolor 3' 'd h.3' 'set ccolor 4' 'd h.4' 'draw title 4 experiments comparison' 'printim compare.png white'

这个脚本把四组试验的 500 百帕高度场叠在一起，颜色区分。h.1到h.4对应四个数据文件里的高度场变量。set lon和set lat限定东北、华北区域，和实习要求的范围一致。printim输出 PNG，方便贴进报告。

参数说明：set gxout contour是等值线图，如果要看风场，改成set gxout vector再d u;v。set ccolor后面的数字是颜色号，1 到 15 对应不同颜色。sdfopen要求数据是 NetCDF 格式，如果 Fortran 输出的是二进制，需要先用grib2ctl或手写.ctl转成 GrADS 能读的格式。

验证方法：跑完四组试验后，重点看两个指标。一是 24 小时预报场和实况场的低压中心位置偏差，实习报告里提到偏北 5 个纬距左右，你可以用 GrADS 的amin和amax函数自动找低压中心经纬度，和实况对比。二是贝加尔湖低压槽的形态，如果槽线位置和走向基本一致，说明模式对影响中国的主要系统预报能力尚可。

从那以后我每次跑正压原始方程模式，都会先把s设成 0.3 跑一组基准，再逐步调参，绝不一上来就上 0.5。时间步长也强制从 300 秒起步，确认稳定后再翻倍。这套习惯帮我省了不少重跑的时间。希望帮到你。

本文还有配套的精品资源，点击获取