简介这份正压原始方程模式实习报告面向大气科学、气象学及数值天气预报方向的本科生与研究生用于课程实习、上机实验与报告撰写参考。报告以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(l1) then do i2,m-1 do j2,n-1 w(i,j)a(i,j)s*(a(i-1,j)a(i1,j)a(i,j-1)a(i,j1)-4*a(i,j))/4.0 end do end do do i2,m-1 do j2,n-1 a(i,j)w(i,j) end do end do return else do i2,m-1 do j2,n-1 w(i,j)a(i,j)s*(a(i-1,j)a(i1,j)a(i,j-1)a(i,j1)-4*a(i,j))/4.0 end do end do do i2,m-1 do j2,n-1 a(i,j)w(i,j) end do end do do i2,m-1 do j2,n-1 w(i,j)a(i,j)(-s)*(a(i-1,j)a(i1,j)a(i,j-1)a(i,j1)-4*a(i,j))/4.0 end do end do do i2,m-1 do j2,n-1 a(i,j)w(i,j) end do end do endif return end逻辑说明l1时只做一次正平滑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) * ∂Φ/∂yv (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 pm-1 qn-1 do i1,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 j2,q ua(i,j)-rm(i,j)*9.8*(za(i,j1)-za(i,j-1))/(2.0*f(i,j)*d) end do end do do j1,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 i2,p va(i,j)rm(i,j)*9.8*(za(i1,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,j1)-za(i,j-1)除以2dva(i,j)用za(i1,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 -fcheckall-O2开优化-fcheckall打开运行时检查数组越界和未初始化变量会直接报错比事后查玄学问题省事得多。编译通过后把初值数据按代码要求的格式写成文本文件一般是m行n列的位势高度场。3.2 试验一正平滑与正逆平滑对比这个试验改的是ssip子程序的l参数。l1只做正平滑l2做正逆平滑。跑两次分别输出 24 小时预报场用 GrADS 叠在一起看差异。# 正平滑 ./model input_l1.txt output_l1.txt # 正逆平滑 ./model input_l2.txt output_l2.txt正平滑会持续耗散短波场变得过于光滑低压中心强度可能被削弱。正逆平滑在滤掉噪声的同时把振幅拉回来一部分预报场更接近实况。但逆平滑不是万能的如果s太大正逆平滑也会引入新的计算噪声。我一般先用s0.5跑一组看预报场和实况的偏差再微调。3.3 试验二地转风差分格式对比这个试验改的是cgw子程序里边界点和内点的差分方式。代码里已经写了前差和后差两种边界处理内点固定中央差分。你可以把边界点的单侧差分改成中央差分需要虚拟点或者把内点的中央差分改成前差/后差观察初值风场的变化。! 把内点中央差分改成前差 ua(i,j)-rm(i,j)*9.8*(za(i,j1)-za(i,j))/(f(i,j)*d) va(i,j)rm(i,j)*9.8*(za(i1,j)-za(i,j))/(f(i,j)*d)前差会让风场整体偏移后差偏移方向相反中央差分最平衡。边界点用单侧差分是必须的因为边界外没有数据。这个试验的价值在于让你直观看到差分格式对初值场的敏感性——初值差一点24 小时后预报场可能差很多。3.4 试验三与试验四边界平滑和时间平滑的开关边界平滑控制的是ssip子程序里循环的范围。默认只平滑内点边界点不动。如果把循环改成do i1,m和do j1,n边界点也参与平滑但固定侧边界条件要求边界值不变所以平滑后要把边界值还原。时间平滑是在时间积分过程中每隔几步做一次平滑抑制高频振荡。代码里没有直接给时间平滑的子程序常见做法是在时间步循环里插入ssip调用每 6 步或 12 步平滑一次。! 时间步循环内插入平滑 if (mod(nt, 6) 0) then call ssip(h, w, 0.5, m, n, 1, 1) end ifnt是时间步计数每 6 步做一次正平滑。不做时间平滑时24 小时积分后可能出现明显的短波噪声高度场等值线变得锯齿状。做了时间平滑场更干净但过度平滑会抹掉真实的中尺度系统。4. 避坑与排查这份代码跑不起来时先看这几条4.1 数组越界导致段错误现象编译通过运行几秒后报Segmentation fault没有任何其他提示。原因ssip和cgw子程序里循环边界写的是2,m-1和2,n-1但如果主程序传入的m,n和数组实际维度不一致或者初值数据行列数读错就会越界。解决在编译时加-fcheckall运行时会打印越界的具体位置。同时检查读初值的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 秒起步确认稳定后再翻倍。这套习惯帮我省了不少重跑的时间。希望帮到你。本文还有配套的精品资源点击获取