简介本资源是一份面向气象、大气科学及相关专业高年级本科生或研究生的数值天气预报实践教学材料聚焦正压原始方程模式的核心原理与编程实现解决理论学习向工程实践转化的关键训练需求。文档完整呈现了以1973年4月29日东北—华北500hPa实测场为初值的24小时有限区域预报全流程涵盖五点平滑子程序含正/逆平滑对比、地转风初值子程序含前差/后差/中心差分格式试验、边界与时间平滑敏感性分析等4组关键数值试验附详细Fortran代码、计算框图、预报场图形及偏差分析。资源为单个Word文档.doc大小361KB结构清晰、注释详实便于代码复现、结果比对与教学讲解。目前已有766人学习下载是掌握数值模式编程、理解守恒平流格式、提升GrADS绘图与Fortran调试能力的典型实习范本。1. 这份 Fortran 实习报告不是“过时文档”而是数值天气预报最硬核的入门切口很多人看到“正压原始方程模式实习报告.doc”第一反应是这不就是上世纪70年代的老古董用 Fortran 写、手绘格点图、连 Python 都没影子——现在谁还这么干但恰恰相反这份报告里藏着现代数值天气预报系统最底层的逻辑骨架。它不讲深度学习或大模型而是用 20 行核心 Fortran 代码把位势高度平流、地转风生成、边界处理、时间步进这些不可绕过的物理约束全部压缩进一个有限区域、500 hPa、24 小时的闭环计算中。你不需要部署 GFS 或 ECMWF只要在本地装好 gfortran GrADS就能复现从初值构造→平滑滤波→时间积分→场量输出的完整链路。它面向的是气象/大气科学专业高年级本科生和刚接触模式开发的研究生——不是教你怎么调参而是逼你亲手推导差分格式、调试数组越界、验证守恒性、比对预报偏差。如果你正在学《数值天气预报》课程或者准备参与 WRF/MPAS 的二次开发这份报告不是历史资料而是你第一个能真正“跑通”的可调试、可修改、可验证的原始方程最小可行实现。2. 正压原始方程的物理约束与 Fortran 实现选择为什么必须用五点平滑和地转风初值2.1 正压假设下的动力学简化从连续方程到离散迭代正压原始方程Barotropic Primitive Equations并非简化版“玩具模型”而是在特定尺度下对大气运动的合理近似。其核心假设是垂直方向密度均匀ρ const因此位势高度 Φ 与气压 P 线性相关Φ g·z ≈ (R/P₀)·T·ln(P₀/P)且水平风场完全由地转平衡主导。此时控制方程退化为仅含位势高度 H单位gpm和水平风分量 u、v 的两个方程连续方程质量守恒∂H/∂t ∇·(HV) 0动量方程地转近似平流项∂u/∂t −u∂u/∂x − v∂u/∂y − f·v (1/H)∂(Hv)/∂y∂v/∂t −u∂v/∂x − v∂v/∂y f·u − (1/H)∂(Hu)/∂x注意这里没有温度方程、没有垂直运动 w、没有非绝热加热项。所有预报变量都只依赖于水平二维格点i,j和时间层 k。这种简化使计算量降低两个数量级但保留了中纬度西风带、槽脊移动、涡度平流等关键动力过程。实习要求中指定“二次守恒平流格式”即采用 Arakawa A-grid 上的二次中心差分如 (u_{i1,j}−u_{i−1,j})/(2Δx)并强制满足离散形式的动能与涡度守恒——这是避免计算崩溃的底线要求而非可选项。提示很多初学者误以为“正压简单”实则正压模式对初值敏感度极高。若初始风场不严格满足地转平衡即 |∇×V| ≈ f后续积分会迅速激发重力波噪声导致位势高度场在几小时后出现高频振荡。这就是实习强制要求“地转风初值子程序”的根本原因——它不是辅助模块而是数值稳定的前置闸门。2.2 五点平滑子程序的算法本质与 Fortran 实现细节五点平滑Five-point smoother在气象数值模式中承担双重角色一是抑制因差分离散引入的 2Δx 尺度虚假振荡即“计算模”二是模拟未解析尺度的湍流耗散效应。其实质是应用一个加权移动平均核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其中s是平滑系数通常取 0.25~0.5括号内为拉普拉斯算子的四邻域离散近似。该公式等价于w(i,j) (1−s)·a(i,j) s/4·[a(i−1,j)a(i1,j)a(i,j−1)a(i,j1)]即中心点权重为(1−s)四个邻点各占s/4。当s0.5时权重分布为[0.5, 0.125, 0.125, 0.125, 0.125]符合经典五点平滑定义。实习报告中的ssip子程序通过l参数控制执行模式l1仅做一次正向平滑forward smoothingl/1先正向平滑再以−s系数做逆向平滑即w(i,j) a(i,j) − s·∇²a构成“正逆平滑”组合这种设计并非随意——正向平滑抑制高频噪声但会模糊锋区逆向平滑则增强梯度类似锐化二者组合可在保特征前提下提升信噪比。Fortran 实现中需特别注意数组索引边界do i2,m−1; do j2,n−1明确排除了第1行/列和末行/列避免a(i−1,j)访问a(0,j)导致段错误。这也是为何报告强调“固定水平侧边界条件”边界格点不参与平滑其值由外部约束如嵌套边界或周期性延拓给定。subroutine ssip(a,w,s,m,n,k,l) implicit none integer m,n,k,l,i,j real a(m,n),w(m,n),s if(l .eq. 1) then ! 正向平滑抑制噪声 do i 2, m-1 do j 2, n-1 w(i,j) a(i,j) s*(a(i-1,j)a(i1,j)a(i,j-1)a(i,j1)-4.0*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 else ! 正逆组合先正向再逆向 do i 2, m-1 do j 2, n-1 w(i,j) a(i,j) s*(a(i-1,j)a(i1,j)a(i,j-1)a(i,j1)-4.0*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(i1,j)a(i,j-1)a(i,j1)-4.0*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 subroutine ssip参数说明a(m,n)输入/输出的二维场如位势高度 Hw(m,n)工作数组用于暂存中间结果s平滑强度系数s0关闭平滑s0.5为常用值m,n格点数x,y方向必须 ≥5 否则循环无效l控制标志l1为单次正向l≠1为正逆组合常见错误排查若编译时报错Segmentation fault (core dumped)90% 源于m,n设置过小如m3导致i2,m−1循环上限为1i2超出范围或a数组未正确分配内存。建议在主程序中添加print *, Grid size: , m, n验证输入。2.3 地转风初值子程序的物理推导与差分格式选择地转风Geostrophic Wind是正压模式初值构建的物理基石。其理论公式为u_g −(g/f)·∂Φ/∂y,v_g (g/f)·∂Φ/∂x其中f 2Ωsinφ为科里奥利参数g9.8 m/s²Φ为位势高度。实习报告中使用的cgw子程序将此公式离散化并针对边界格点采用不同差分策略格点位置u 分量计算方式v 分量计算方式物理含义左右边界j1, jn单侧后差/前差中心差i2..m−1避免越界牺牲精度保存在上下边界i1, im中心差j2..n−1单侧后差/前差同上内部格点i2..m−1, j2..n−1二阶中心差二阶中心差最高精度具体实现中ua(i,1)和ua(i,n)使用一阶后差/前差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)而内部点ua(i,j)使用二阶中心差ua(i,j) −rm(i,j)*9.8*(za(i,j1)−za(i,j−1))/(2.0*f(i,j)*d)同理va在 x 方向边界用单侧差分内部用中心差。这里的d是格距单位米rm是平均空气密度单位kg/m³f是科氏参数单位s⁻¹。所有变量均为real类型implicit none强制显式声明杜绝隐式类型错误。注意cgw子程序中pm−1,qn−1的设定是为了在do i2,p循环中自然覆盖i2到im−1避免im越界访问za(i1,j)。这是 Fortran 数值编程的经典边界处理技巧比直接写do i2,m−1更易维护。3. 四组数值试验的设计逻辑与可复现实操步骤从对比到归因3.1 正平滑 vs 正逆平滑如何量化平滑策略对预报误差的影响这两组试验直指数值稳定性与物理保真度的权衡。正平滑l1仅做一次低通滤波会平抑所有小尺度扰动包括真实的锋面梯度正逆平滑l≠1先平滑再“反平滑”相当于对拉普拉斯算子做两次操作∇²(∇²a)其频谱响应在中尺度有轻微增强有利于维持槽脊结构。要复现该对比需在主程序中调用ssip两次# 编译并运行两种配置假设主程序为 main.f gfortran -o forecast_forward main.f ssip.f cgw.f gfortran -o forecast_forward_inverse main.f ssip.f cgw.f关键修改在主程序调用处! 正平滑配置试验① call ssip(H, W, 0.35, m, n, 1, 1) ! l1 ! 正逆平滑配置试验② call ssip(H, W, 0.35, m, n, 1, 2) ! l2预报结果验证方法空间误差计算预报场与实况场教材图的均方根误差RMSERMSE sqrt( sum[(H_f(i,j)−H_obs(i,j))²] / (m×n) )系统性偏差统计低压中心位置偏移量纬距/经距结构保真度提取沿 45°N 的位势高度剖面对比槽深gpm和槽宽格点数实习报告图7/8显示正逆平滑的低压中心东移距离更接近实况偏移约3个纬距 vs 正平滑的5个纬距证明其对涡度平流的表征更优。这并非偶然——正逆组合实质上逼近了双调和滤波biharmonic filter对中尺度系统能量耗散更少。3.2 地转风差分格式试验前差、后差、中心差的精度代价分析该试验暴露初值质量对模式性能的决定性影响。三种差分格式的截断误差阶数不同前差/后差O(Δx)一阶精度边界适用但相位误差大中心差O(Δx²)二阶精度内部最优但边界需特殊处理在cgw.f中可通过修改ua和va的边界计算式来切换! 替换 ua(i,1) 的后差为前差仅用于试验 ua(i,1) -rm(i,1)*9.8*(za(i,3)-za(i,1))/(2.0*f(i,1)*d) ! 伪中心差需 za(i,3)但实际操作中za(i,3)可能不存在故试验②本质是验证当被迫使用低阶差分时模式是否仍能维持24小时预报可用性结果图8/9表明前差初值导致预报场在6小时后即出现虚假西风急流而后差初值虽略偏弱但结构更稳定。这印证了数值分析结论初值误差会随时间指数放大Lyapunov 指数因此高阶差分不仅是精度问题更是稳定性门槛。3.3 边界平滑与时间平滑试验识别模式“失稳源”的诊断工具这两个试验针对两类典型失稳机制边界平滑关闭试验③若不平滑边界格点模式在侧边界处易激发出反射重力波表现为位势高度场在边界附近出现同心圆状振荡见图10。解决方案是在ssip调用前对i1,m和j1,n行列单独做1D平滑或采用海绵边界sponge layer。时间平滑关闭试验④时间平滑如 Robert-Asselin 滤波用于抑制时间积分中的计算模。关闭后H场在第3-5预报小时会出现高频“抖动”振幅达50 gpm远超真实天气变率。此时需检查时间步长Δt是否满足 CFL 条件Δt Δx / max(|u|,|v|)。实习中Δt通常设为 300 秒5分钟若风速达 20 m/sΔx至少需 6 km。验证命令Linux 下# 提取第12小时预报场的边界行i1 grads EOF open forecast.ctl set t 12 set z 1 set x 1 100 set y 1 1 d hgt quit EOF # 输出为 ASCII用 awk 统计标准差 awk {sum\$1; sumsq\$1*\$1} END {print STD:, sqrt(sumsq/NR - (sum/NR)^2)} grads.dat若 STD 10 gpm即判定边界失稳。4. GrADS 可视化与误差归因从图形对比到物理机制诊断4.1 GrADS 控制文件.ctl编写规范与常见陷阱GrADS 是本实习唯一指定绘图工具其.ctl文件定义了数据的时空结构。一个典型forecast.ctl应包含DSET ^forecast.dat TITLE 500hPa Height and Wind Forecast UNDEF -9999 XDEF 100 LINEAR 110 0.5 # 100格点起始经度110E格距0.5° YDEF 80 LINEAR 30 0.5 # 80格点起始纬度30N格距0.5° ZDEF 1 LEVELS 500 TDEF 49 LINEAR 29:08Z04jan1973 1HR # 49个时次0h,1h,...,48h VARS 3 hgt 0 0 Z500 500hPa geopotential height (gpm) u 0 0 Z500 500hPa u-wind (m/s) v 0 0 Z500 500hPa v-wind (m/s) ENDVARS关键陷阱DSET路径必须为相对路径且forecast.dat必须是二进制大端序Big-endian格式。Fortran 默认写入大端序但若在 x86_64 Linux 编译需加-fconvertbig-endiangfortran -fconvertbig-endian -o main main.fTDEF的起始时间必须与初值时间严格一致1973年4月29日08时否则set t 24会定位到错误时刻。VARS中变量顺序必须与 Fortranwrite语句顺序完全匹配否则d hgt会读取u数据。4.2 误差归因的三层诊断法从图形到方程实习报告图12指出“低压中心偏北5个纬距”这不能止步于现象描述而应归因到物理过程。推荐按以下三层递进分析第一层场量诊断GrADS 直接计算# 计算涡度 advection-u*∂ζ/∂x - v*∂ζ/∂y set x 1 100; set y 1 80; set t 24 define zeta (v[i1,j]-v[i-1,j] - u[i,j1]u[i,j-1])/(2*0.5*111) # 简化涡度 define adv -u*d(zeta,x) - v*d(zeta,y) d adv若预报场中adv在低压中心西侧为负冷平流而实况为正暖平流则说明模式低估了槽前西南气流的输送。第二层方程残差诊断修改 Fortran 源码在main.f的时间循环中插入残差计算! 在每个时间步后计算连续方程残差 resid 0.0 do i2,m-1 do j2,n-1 dHdt (Hnew(i,j)-Hold(i,j))/dt div (H(i1,j)*u(i1,j)-H(i-1,j)*u(i-1,j) H(i,j1)*v(i,j1)-H(i,j-1)*v(i,j-1))/(2*dx) resid resid abs(dHdt div) end do end do print *, Time step, k, Residual:, resid/(m*n)若resid在第10步后突增10倍说明平流项离散格式在该区域失效。第三层初值敏感性试验扰动法对初值za添加随机扰动δza 0.1 * rand()重新运行48小时。若预报场 RMSE 增幅 30%则证实初值误差是主要误差源——这正是正压模式的固有局限也解释了为何现代业务模式必须耦合资料同化。提示当 GrADS 报错fortran显示无法启动程序95% 源于.dat文件字节长度不匹配。用ls -l forecast.dat查看大小应等于m*n*4*3*493变量×49时次×每变量4字节。若不符检查 Fortranwrite是否遗漏rec参数或format错误。5. Fortran 模式调试的黄金三原则从段错误到物理一致性5.1 编译期防御用 gfortran 的静态检查堵住 70% 的错误现代 gfortran 提供强大静态分析能力应在编译阶段启用gfortran -Wall -Wextra -fcheckall -fbounds-check \ -finit-realnan -finit-integer-999 \ -o debug_mode ssip.f cgw.f main.f参数详解-Wall -Wextra开启所有警告如未初始化变量、无用赋值-fcheckall运行时检查数组越界、指针空引用、递归调用-fbounds-check强制检查a(i,j)的i,j是否在声明范围内-finit-realnan将未初始化real变量设为 NaN避免静默错误-finit-integer-999同理初始化整型变量若程序运行时报Program received signal SIGFPE: Floating-point exception立即检查f(i,j)是否为0赤道附近f0会导致除零或d是否被误设为0。5.2 运行时验证三个必查的物理守恒量正压模式必须满足三项基本守恒可在每个时间步后插入验证! 1. 总位势高度守恒质量守恒 sum_H 0.0 do i1,m; do j1,n; sum_H sum_H H(i,j); end do; end do print *, Step, k, Total H:, sum_H ! 2. 总动能守恒无外力时 sum_KE 0.0 do i1,m; do j1,n; sum_KE sum_KE 0.5*(u(i,j)**2 v(i,j)**2); end do; end do print *, Step, k, Total KE:, sum_KE ! 3. 总绝对涡度守恒无摩擦时 sum_zeta 0.0 do i2,m-1; do j2,n-1 zeta (v(i1,j)-v(i-1,j) - u(i,j1)u(i,j-1))/(2*dx) sum_zeta sum_zeta zeta end do; end do print *, Step, k, Total Zeta:, sum_zeta若sum_H在24小时内变化 0.1%说明平流格式未严格守恒需检查ssip是否误改了边界格点。5.3 图形级调试用 GrADS 的query命令定位异常格点当预报场出现局部“炸点”如某格点hgt9999不要盲目重跑用 GrADS 快速定位open forecast.ctl set t 12 set z 1 query dims # 查看当前维度范围 query file # 查看变量信息 d hgt # 显示全场 q gxvals # 输出当前显示区域所有格点值若q gxvals输出中某行含9999.000记下其i,j坐标如i45,j32然后回到 Fortran 源码搜索H(45,32)的所有赋值语句——大概率是cgw中f(45,32)为0或d0导致除零或ssip循环未覆盖该点。最后提醒这份实习报告的价值不在于它“古老”而在于它剔除了所有现代模式的工程封装MPI 并行、NetCDF I/O、物理参数化让你直面原始方程的数学内核。当你亲手修复一个Array bound violation错误或发现s0.4比s0.3更好地平衡了噪声与分辨率你就真正跨过了数值天气预报的第一道门槛——不是调包而是造轮子。本文还有配套的精品资源点击获取