简介：本资源是一份面向气象、大气科学及相关专业高年级本科生或研究生的数值天气预报实践教学材料，聚焦正压原始方程模式的核心原理与编程实现，解决理论学习向工程实践转化的关键训练需求。文档完整呈现了以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_{i+1,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(i+1,j)+a(i,j−1)+a(i,j+1) − 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(i+1,j)+a(i,j−1)+a(i,j+1)]
即中心点权重为(1−s)，四个邻点各占s/4。当s=0.5时，权重分布为[0.5, 0.125, 0.125, 0.125, 0.125]，符合经典五点平滑定义。
实习报告中的ssip子程序通过l参数控制执行模式：
l==1：仅做一次正向平滑（forward smoothing）l/=1：先正向平滑，再以−s系数做逆向平滑（即w(i,j) = a(i,j) − s·∇²a），构成“正逆平滑”组合
这种设计并非随意——正向平滑抑制高频噪声但会模糊锋区，逆向平滑则增强梯度（类似锐化），二者组合可在保特征前提下提升信噪比。Fortran 实现中需特别注意数组索引边界：do i=2,m−1; do j=2,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(i+1,j)+a(i,j-1)+a(i,j+1)-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(i+1,j)+a(i,j-1)+a(i,j+1)-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(i+1,j)+a(i,j-1)+a(i,j+1)-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)：输入/输出的二维场（如位势高度 H）w(m,n)：工作数组，用于暂存中间结果s：平滑强度系数，s=0关闭平滑，s=0.5为常用值m,n：格点数（x,y方向），必须 ≥5 否则循环无效l：控制标志，l=1为单次正向，l≠1为正逆组合
常见错误排查：若编译时报错Segmentation fault (core dumped)，90% 源于m,n设置过小（如m=3导致i=2,m−1循环上限为1，i=2超出范围）或a数组未正确分配内存。建议在主程序中添加print *, 'Grid size: ', m, n验证输入。
2.3 地转风初值子程序的物理推导与差分格式选择
地转风（Geostrophic Wind）是正压模式初值构建的物理基石。其理论公式为：u_g = −(g/f)·∂Φ/∂y,v_g = (g/f)·∂Φ/∂x
其中f = 2Ωsinφ为科里奥利参数，g=9.8 m/s²，Φ为位势高度。实习报告中使用的cgw子程序将此公式离散化，并针对边界格点采用不同差分策略：
| 格点位置 | u 分量计算方式 | v 分量计算方式 | 物理含义 |
|---|---|---|---|
| 左右边界（j=1, j=n） | 单侧后差/前差 | 中心差（i=2..m−1） | 避免越界，牺牲精度保存在 |
| 上下边界（i=1, i=m） | 中心差（j=2..n−1） | 单侧后差/前差 | 同上 |
| 内部格点（i=2..m−1, j=2..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,j+1)−za(i,j−1))/(2.0*f(i,j)*d)
同理va在 x 方向边界用单侧差分，内部用中心差。这里的d是格距（单位：米），rm是平均空气密度（单位：kg/m³），f是科氏参数（单位：s⁻¹）。所有变量均为real类型，implicit none强制显式声明，杜绝隐式类型错误。
注意：
cgw子程序中p=m−1,q=n−1的设定，是为了在do i=2,p循环中自然覆盖i=2到i=m−1，避免i=m越界访问za(i+1,j)。这是 Fortran 数值编程的经典边界处理技巧，比直接写do i=2,m−1更易维护。
3. 四组数值试验的设计逻辑与可复现实操步骤：从对比到归因
3.1 正平滑 vs 正逆平滑：如何量化平滑策略对预报误差的影响？
这两组试验直指数值稳定性与物理保真度的权衡。正平滑（l=1）仅做一次低通滤波，会平抑所有小尺度扰动，包括真实的锋面梯度；正逆平滑（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) ! l=1 ! 正逆平滑配置（试验②） call ssip(H, W, 0.35, m, n, 1, 2) ! l=2预报结果验证方法：
- 空间误差：计算预报场与实况场（教材图）的均方根误差（RMSE）
RMSE = 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调用前，对i=1,m和j=1,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小时预报场的边界行（i=1） 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 Z=500 500hPa geopotential height (gpm) u 0 0 Z=500 500hPa u-wind (m/s) v 0 0 Z=500 500hPa v-wind (m/s) ENDVARS关键陷阱：
DSET路径必须为相对路径，且forecast.dat必须是二进制大端序（Big-endian）格式。Fortran 默认写入大端序，但若在 x86_64 Linux 编译，需加-fconvert=big-endian：gfortran -fconvert=big-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[i+1,j]-v[i-1,j] - u[i,j+1]+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 i=2,m-1 do j=2,n-1 dHdt = (Hnew(i,j)-Hold(i,j))/dt div = (H(i+1,j)*u(i+1,j)-H(i-1,j)*u(i-1,j) + & H(i,j+1)*v(i,j+1)-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*49（3变量×49时次×每变量4字节）。若不符，检查 Fortranwrite是否遗漏rec=参数或format错误。
5. Fortran 模式调试的黄金三原则：从段错误到物理一致性
5.1 编译期防御：用 gfortran 的静态检查堵住 70% 的错误
现代 gfortran 提供强大静态分析能力，应在编译阶段启用：
gfortran -Wall -Wextra -fcheck=all -fbounds-check \ -finit-real=nan -finit-integer=-999 \ -o debug_mode ssip.f cgw.f main.f参数详解：
-Wall -Wextra：开启所有警告，如未初始化变量、无用赋值-fcheck=all：运行时检查数组越界、指针空引用、递归调用-fbounds-check：强制检查a(i,j)的i,j是否在声明范围内-finit-real=nan：将未初始化real变量设为 NaN，避免静默错误-finit-integer=-999：同理初始化整型变量
若程序运行时报Program received signal SIGFPE: Floating-point exception，立即检查f(i,j)是否为0（赤道附近f=0会导致除零），或d是否被误设为0。
5.2 运行时验证：三个必查的物理守恒量
正压模式必须满足三项基本守恒，可在每个时间步后插入验证：
! 1. 总位势高度守恒（质量守恒） sum_H = 0.0 do i=1,m; do j=1,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 i=1,m; do j=1,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 i=2,m-1; do j=2,n-1 zeta = (v(i+1,j)-v(i-1,j) - u(i,j+1)+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命令定位异常格点
当预报场出现局部“炸点”（如某格点hgt=9999），不要盲目重跑，用 GrADS 快速定位：
open forecast.ctl set t 12 set z 1 query dims # 查看当前维度范围 query file # 查看变量信息 d hgt # 显示全场 q gxvals # 输出当前显示区域所有格点值若q gxvals输出中某行含9999.000，记下其i,j坐标（如i=45,j=32），然后回到 Fortran 源码搜索H(45,32)的所有赋值语句——大概率是cgw中f(45,32)为0或d=0导致除零，或ssip循环未覆盖该点。
最后提醒：这份实习报告的价值，不在于它“古老”，而在于它剔除了所有现代模式的工程封装（MPI 并行、NetCDF I/O、物理参数化），让你直面原始方程的数学内核。当你亲手修复一个Array bound violation错误，或发现s=0.4比s=0.3更好地平衡了噪声与分辨率，你就真正跨过了数值天气预报的第一道门槛——不是调包，而是造轮子。
本文还有配套的精品资源，点击获取
