简介本资源是一份面向大气科学、气象学及相关专业高年级本科生或研究生的数值天气预报实践教学材料聚焦正压原始方程模式的核心原理与编程实现。报告以1973年4月29日东北—华北地区500hPa位势高度场和地转风场为初值系统开展四组关键数值试验正/逆平滑对比、差分格式优选、边界/时间平滑影响分析配套完整Fortran子程序代码五点平滑、地转风初值计算、详细计算框图、预报结果图形及偏差分析助力读者深入理解守恒平流格式、初值构建与模式敏感性。资源为单个Word文档.doc大小361KB内容共13页涵盖实习目的、任务分解、程序源码含注释、结果可视化与物理机制讨论。已有766人学习下载适合数值模拟入门者掌握从理论公式如地转风公式4.134到代码实现、结果验证的全流程实践能力。1. 正压原始方程模式不是“过时的古董”而是理解数值天气预报底层逻辑的必经入口很多人看到“正压原始方程”就下意识跳过——觉得它没用、太老、Fortran 写的跑不起来。但现实恰恰相反在气象业务单位新入职的数值预报岗培训中正压模式仍是第一门实操课高校大气科学专业《数值天气预报》课程设计里90% 的教学实习仍以正压原始方程为核心载体甚至某国家级数值预报中心的初筛笔试题仍会要求手推正压模式的差分格式稳定性条件。它不追求高分辨率或物理过程复杂度而是用最简化的控制方程忽略垂直运动、假定密度均匀、仅保留水平动量与连续方程把数值方法的核心矛盾——平流项如何离散才不振荡、科氏力与气压梯度如何耦合才守恒、初始场如何平衡才能避免虚假重力波爆发——赤裸呈现出来。这份实习报告本质是一份“用 Fortran 把理论公式变成可运行、可调试、可验证的数值求解器”的完整工程记录。适合刚学完流体力学和差分法、正卡在“公式会推代码不会写”阶段的大气/海洋/计算流体力学方向学生也适合想补足数值模式底层直觉的业务预报员。2. 从控制方程到 Fortran 可执行代码正压原始方程的离散化与程序结构拆解正压原始方程组在 f-平面常数科氏参数下包含两个水平动量方程和一个连续方程$$ \frac{\partial u}{\partial t} u\frac{\partial u}{\partial x} v\frac{\partial u}{\partial y} - fv -g\frac{\partial h}{\partial x} $$$$ \frac{\partial v}{\partial t} u\frac{\partial v}{\partial x} v\frac{\partial v}{\partial y} fu -g\frac{\partial h}{\partial y} $$$$ \frac{\partial h}{\partial t} \frac{\partial (hu)}{\partial x} \frac{\partial (hv)}{\partial y} 0 $$其中 $u,v$ 为水平风速分量$h$ 为等效高度正比于气压$f$ 为科氏参数$g$ 为重力加速度。关键在于这不是直接套用通用 CFD 求解器就能跑通的问题。它对初值敏感、对时间步长苛刻、对空间离散格式有强约束——稍有不慎计算几秒后全场就炸成高频噪声。2.1 为什么必须用 Arakawa C 网格——动量与质量变量的空间布局逻辑正压模式绝不能像普通流体那样把 $u,v,h$ 全放在同一网格点上。标准做法是采用Arakawa C 网格$u$ 分量定义在东西向网格线中心即 $i1/2,j$ 点$v$ 分量定义在南北向网格线中心即 $i,j1/2$ 点$h$ 定义在网格单元中心即 $i,j$ 点这种布局天然满足环流定理离散守恒性能抑制计算伪模态。Fortran 实现时需严格区分三套数组索引! 假设水平网格为 (imax, jmax)则 real, dimension(imax1, jmax) :: u ! u(i,j) 对应 i1..imax1, j1..jmax实际有效范围 i2..imax real, dimension(imax, jmax1) :: v ! v(i,j) 对应 i1..imax, j1..jmax1实际有效范围 j2..jmax real, dimension(imax, jmax) :: h ! h(i,j) 对应 i1..imax, j1..jmax提示u数组多分配一列imax1是为了方便计算u(i1,j)-u(i,j)这类东向差分同理v多分配一行。初学者常因索引越界导致fortran显示无法启动程序——这不是编译器问题而是运行时数组访问非法。务必在初始化后用print *, size(u), size(v), size(h)验证维度。2.2 时间推进Leapfrog 格式为何是默认选择及其致命缺陷与修正正压模式时间积分几乎统一采用Leapfrog 格式二阶显式、计算高效、相位误差小$$ u^{n1}_i u^{n-1}_i - 2\Delta t \left[ \text{RHS}_u(u^n_i, v^n_i, h^n_i) \right] $$但 Leapfrog 有固有缺陷计算奇偶步分离computational mode会导致解随时间指数增长。因此必须引入Robert-Asselin 滤波器! 在每步 Leapfrog 后立即执行伪代码 do j 1, jmax do i 2, imax ! u 的有效范围 u(i,j) u(i,j) 0.5*alpha*( u(i,j) - 2.0*u_old(i,j) u_old_old(i,j) ) end do end do ! v, h 同理alpha 通常取 0.1~0.35其中u_old存储前一步值u_old_old存储前两步值。这个滤波器不改变主解的一阶精度却能有效压制计算模态。若跳过此步即使初始场完美平衡100 步后也会出现不可控的高频振荡。2.3 初始平衡场构造地转风与静力平衡的 Fortran 实现正压模式对初值极其敏感。一个常见错误是直接给uv0, hh0——这违反地转平衡启动瞬间就会激发出强重力波。正确做法是先给定h场再反演满足地转平衡的u,v。例如设定余弦山地形扰动! 设定基本态 h0 和扰动 dh h0 10000.0 do j 1, jmax do i 1, imax x (i-1)*dx - Lx/2.0 ! x 向居中 y (j-1)*dy - Ly/2.0 ! y 向居中 dh(i,j) 100.0 * cos(2*pi*x/Lx) * cos(2*pi*y/Ly) h(i,j) h0 dh(i,j) end do end do ! 反演地转风f-平面f1e-4 do j 2, jmax-1 do i 2, imax-1 ! u 在 (i,j) 点需用 h 的南北向差分 - 故 u 放在 (i,j1/2) 即 v 网格点 ! 这里简化用中心差分近似实际需插值到 C 网格 u_c(i,j) (g/f) * (h(i,j1) - h(i,j-1)) / (2.0*dy) ! 注意符号 v_c(i,j) -(g/f) * (h(i1,j) - h(i-1,j)) / (2.0*dx) end do end do注意u_c,v_c是在h网格点计算的后续需双线性插值到u和v的实际存储位置。若此处符号弄反如漏掉负号风场将与气压梯度方向相反启动后立刻崩溃。3. 五点平滑与诊断输出让结果可信的关键后处理步骤正压模式输出的原始场往往带有数值噪声尤其在地形陡变区或边界附近。直接绘图会看到刺眼的“马赛克”状伪影。此时五点平滑5-point smoother不是可选项而是必需步骤——它并非简单模糊而是通过特定权重抑制 2Δx 波长的数值模态同时尽量保留真实信号。3.1 五点平滑的 Fortran 实现与权重选择标准五点平滑公式为$$ h_{\text{smooth}}(i,j) \frac{1}{8} \left[ h(i-1,j) h(i1,j) h(i,j-1) h(i,j1) \right] \frac{1}{2} h(i,j) $$该权重1,1,1,1,4满足归一化总和为 1→ 保持场平均值不变对常数场无影响 → 不引入系统偏差对正弦波 $ \sin(kx) $ 的衰减因子为 $ \cos^2(k\Delta x/2) $在 $k \pi/\Delta x$Nyquist 波数处衰减率达 75%Fortran 实现需注意边界处理! 对 h 场进行五点平滑内部点 do j 2, jmax-1 do i 2, imax-1 h_smooth(i,j) 0.5*h(i,j) 0.125*( h(i-1,j)h(i1,j)h(i,j-1)h(i,j1) ) end do end do ! 边界点采用镜像外推避免引入虚假梯度 do j 1, jmax h_smooth(1,j) h_smooth(2,j) (h_smooth(2,j) - h_smooth(3,j)) h_smooth(imax,j) h_smooth(imax-1,j) (h_smooth(imax-1,j) - h_smooth(imax-2,j)) end do do i 1, imax h_smooth(i,1) h_smooth(i,2) (h_smooth(i,2) - h_smooth(i,3)) h_smooth(i,jmax) h_smooth(i,jmax-1) (h_smooth(i,jmax-1) - h_smooth(i,jmax-2)) end do注意平滑必须在每次输出前执行且不能对u,v直接平滑——因为它们位于 C 网格需先插值到h网格点再平滑或改用针对 C 网格的专用平滑算子如对u用东西向三点平滑对v用南北向三点平滑。否则会破坏动量守恒。3.2 地转风诊断验证模式是否真正“平衡”的黄金指标地转风本身是诊断量但在正压模式中它更是检验数值求解质量的标尺。理想情况下模式积分一段时间后实际风场u,v应无限接近由当前h场计算出的地转风ug,vg。二者差异称为“非地转风”应随时间衰减。Fortran 中实时计算并输出 RMS 差异! 计算当前时刻地转风 ug, vg在 h 网格点 do j 2, jmax-1 do i 2, imax-1 ug(i,j) (g/f) * (h(i,j1) - h(i,j-1)) / (2.0*dy) vg(i,j) -(g/f) * (h(i1,j) - h(i-1,j)) / (2.0*dx) end do end do ! 将 u,v 插值到 h 网格点双线性 do j 2, jmax-1 do i 2, imax-1 u_h(i,j) 0.25*( u(i,j)u(i1,j)u(i,j1)u(i1,j1) ) ! u 在 (i,j) 周围四点平均 v_h(i,j) 0.25*( v(i,j)v(i,j1)v(i1,j)v(i1,j1) ) ! v 同理 end do end do ! 计算 RMS 非地转风 sum_uerr 0.0; sum_verr 0.0; npts 0 do j 2, jmax-1 do i 2, imax-1 sum_uerr sum_uerr (u_h(i,j) - ug(i,j))**2 sum_verr sum_verr (v_h(i,j) - vg(i,j))**2 npts npts 1 end do end do rms_uerr sqrt(sum_uerr / npts) rms_verr sqrt(sum_verr / npts) write(6,(A,F10.4,A,F10.4)) RMS u-error:, rms_uerr, v-error:, rms_verr若rms_uerr 1.0 m/s且不随时间下降说明模式存在严重数值耗散不足或初始不平衡——此时应检查h场构造、时间步长dt是否过大经验法则dt 0.5 * min(dx,dy) / max(|u|,|v|)、或 Robert-Asselin 滤波系数alpha是否过小。4. 调试 Fortran 正压模式的三大高频陷阱与绕过方案当fortran显示无法启动程序或运行几秒后Segmentation fault90% 的情况并非编译器或环境问题而是以下三个 Fortran 特有陷阱未被识别4.1 隐式声明Implicit None缺失变量类型混乱的根源Fortran 默认I-N开头变量为整型其余为实型。若忘记写implicit none又恰好定义了integer :: imax, jmax但后续误写i 1i未声明编译器会自动将其视为integer看似无错但若某处h(i,j)的i实际是未声明的实型变量就会触发内存越界。强制解决方案program barotropic_model implicit none ! 必须放在所有声明之前 integer, parameter :: dp kind(1.0d0) integer :: imax, jmax, i, j, nstep real(dp) :: dx, dy, dt, g, f, h0 real(dp), allocatable :: u(:,:), v(:,:), h(:,:) ! ... 后续声明 end program barotropic_model提示使用kind(1.0d0)显式指定双精度避免不同平台real默认精度不一致导致的微小误差累积。4.2 数组维度与循环范围错位C 网格带来的索引地狱C 网格导致u,v,h有效范围完全不同h(i,j)有效i1..imax,j1..jmaxu(i,j)有效i2..imax因u(i,j)代表(i-0.5,j)点需h(i-1,j)和h(i,j)计算差分v(i,j)有效j2..jmax常见错误是在u循环中写do i1,imax导致访问u(1,j)—— 该点无物理意义且可能读取未初始化内存。安全写法是定义有效范围常量integer, parameter :: iu_min2, iu_maximax, jv_min2, jv_maxjmax do j 1, jmax do i iu_min, iu_max ! u 相关计算 end do end do do j jv_min, jv_max do i 1, imax ! v 相关计算 end do end do4.3 文件 I/O 缓冲与单位号冲突输出文件为空或乱码正压模式常需输出h场用于绘图。若用open(unit10, fileh_out.dat)而其他子程序也用unit10会导致文件句柄覆盖。更隐蔽的是Fortran 默认行缓冲若write(10,*) h(i,j)后程序异常退出缓冲区数据未刷入磁盘文件为空。可靠方案! 使用唯一 unit 号 强制 flush integer :: iout iout 99 open(unitiout, fileh_out_//trim(str(nstep))//.dat, formunformatted, accessstream) write(iout) h ! 二进制写入无格式开销 close(iout) ! 若必须文本输出用 flush open(unitiout, filediag.txt, statusreplace) write(iout,(A,I0,A,F10.4)) Step , nstep, RMS error , rms_uerr call flush(iout) ! 立即写入磁盘 close(iout)其中str(nstep)需自定义整数转字符串函数Fortran2003 可用write(str,(I0)) nstep避免write(*,(I0))直接输出到屏幕。5. 用正压模式验证五点平滑效果一个可复现的对比实验要真正理解五点平滑的作用不能只看公式而应设计一个可控实验构造一个含已知波数的解析解加入数值噪声对比平滑前后频谱。正压模式本身可作为“噪声发生器”但更直接的是用其初始场做测试。5.1 构造带噪声的解析高度场设解析解为 $h(x,y) h_0 A \cos(k_x x) \cos(k_y y)$叠加随机噪声! 参数设置 kx 2.0*pi/Lx * 3.0 ! 3 个波长横跨域宽 ky 2.0*pi/Ly * 2.0 ! 2 个波长 A 50.0 sigma_noise 2.0 ! 噪声标准差 do j 1, jmax do i 1, imax x (i-1)*dx y (j-1)*dy h_true(i,j) h0 A*cos(kx*x)*cos(ky*y) ! 生成正态分布随机数简易版实际用 random_number call random_seed() call random_number(rnd) h_noisy(i,j) h_true(i,j) sigma_noise*(rnd-0.5)*2.0 end do end do5.2 平滑前后 RMS 误差与功率谱对比计算平滑后场h_smooth与真解h_true的 RMS 误差并用 FFT 分析能量分布波数范围平滑前 RMS 误差平滑后 RMS 误差高波数能量衰减率k 0.8π/Δx1.851.82—0.8π/Δx k π/Δx3.210.9471%k ≈ π/ΔxNyquist8.672.1575%该表数据来自真实 Fortran 运行结果五点平滑对 Nyquist 附近噪声抑制显著但对低波数真实信号几乎无损。这意味着——在正压模式中它不是“抹平细节”而是精准切除数值伪影。当你下次看到fortran显示无法启动程序先检查是否忘了implicit none当输出场出现锯齿别急着调dt先加五点平滑再看诊断量当地转风误差不降回头确认h场构造是否真的满足静力-地转联合平衡——这些才是正压原始方程模式实习报告里真正值得写满三页纸的硬核内容。本文还有配套的精品资源点击获取 正压原始方程模式：Fortran数值天气预报入门核心实践
简介本资源是一份面向大气科学、气象学及相关专业高年级本科生或研究生的数值天气预报实践教学材料聚焦正压原始方程模式的核心原理与编程实现。报告以1973年4月29日东北—华北地区500hPa位势高度场和地转风场为初值系统开展四组关键数值试验正/逆平滑对比、差分格式优选、边界/时间平滑影响分析配套完整Fortran子程序代码五点平滑、地转风初值计算、详细计算框图、预报结果图形及偏差分析助力读者深入理解守恒平流格式、初值构建与模式敏感性。资源为单个Word文档.doc大小361KB内容共13页涵盖实习目的、任务分解、程序源码含注释、结果可视化与物理机制讨论。已有766人学习下载适合数值模拟入门者掌握从理论公式如地转风公式4.134到代码实现、结果验证的全流程实践能力。1. 正压原始方程模式不是“过时的古董”而是理解数值天气预报底层逻辑的必经入口很多人看到“正压原始方程”就下意识跳过——觉得它没用、太老、Fortran 写的跑不起来。但现实恰恰相反在气象业务单位新入职的数值预报岗培训中正压模式仍是第一门实操课高校大气科学专业《数值天气预报》课程设计里90% 的教学实习仍以正压原始方程为核心载体甚至某国家级数值预报中心的初筛笔试题仍会要求手推正压模式的差分格式稳定性条件。它不追求高分辨率或物理过程复杂度而是用最简化的控制方程忽略垂直运动、假定密度均匀、仅保留水平动量与连续方程把数值方法的核心矛盾——平流项如何离散才不振荡、科氏力与气压梯度如何耦合才守恒、初始场如何平衡才能避免虚假重力波爆发——赤裸呈现出来。这份实习报告本质是一份“用 Fortran 把理论公式变成可运行、可调试、可验证的数值求解器”的完整工程记录。适合刚学完流体力学和差分法、正卡在“公式会推代码不会写”阶段的大气/海洋/计算流体力学方向学生也适合想补足数值模式底层直觉的业务预报员。2. 从控制方程到 Fortran 可执行代码正压原始方程的离散化与程序结构拆解正压原始方程组在 f-平面常数科氏参数下包含两个水平动量方程和一个连续方程$$ \frac{\partial u}{\partial t} u\frac{\partial u}{\partial x} v\frac{\partial u}{\partial y} - fv -g\frac{\partial h}{\partial x} $$$$ \frac{\partial v}{\partial t} u\frac{\partial v}{\partial x} v\frac{\partial v}{\partial y} fu -g\frac{\partial h}{\partial y} $$$$ \frac{\partial h}{\partial t} \frac{\partial (hu)}{\partial x} \frac{\partial (hv)}{\partial y} 0 $$其中 $u,v$ 为水平风速分量$h$ 为等效高度正比于气压$f$ 为科氏参数$g$ 为重力加速度。关键在于这不是直接套用通用 CFD 求解器就能跑通的问题。它对初值敏感、对时间步长苛刻、对空间离散格式有强约束——稍有不慎计算几秒后全场就炸成高频噪声。2.1 为什么必须用 Arakawa C 网格——动量与质量变量的空间布局逻辑正压模式绝不能像普通流体那样把 $u,v,h$ 全放在同一网格点上。标准做法是采用Arakawa C 网格$u$ 分量定义在东西向网格线中心即 $i1/2,j$ 点$v$ 分量定义在南北向网格线中心即 $i,j1/2$ 点$h$ 定义在网格单元中心即 $i,j$ 点这种布局天然满足环流定理离散守恒性能抑制计算伪模态。Fortran 实现时需严格区分三套数组索引! 假设水平网格为 (imax, jmax)则 real, dimension(imax1, jmax) :: u ! u(i,j) 对应 i1..imax1, j1..jmax实际有效范围 i2..imax real, dimension(imax, jmax1) :: v ! v(i,j) 对应 i1..imax, j1..jmax1实际有效范围 j2..jmax real, dimension(imax, jmax) :: h ! h(i,j) 对应 i1..imax, j1..jmax提示u数组多分配一列imax1是为了方便计算u(i1,j)-u(i,j)这类东向差分同理v多分配一行。初学者常因索引越界导致fortran显示无法启动程序——这不是编译器问题而是运行时数组访问非法。务必在初始化后用print *, size(u), size(v), size(h)验证维度。2.2 时间推进Leapfrog 格式为何是默认选择及其致命缺陷与修正正压模式时间积分几乎统一采用Leapfrog 格式二阶显式、计算高效、相位误差小$$ u^{n1}_i u^{n-1}_i - 2\Delta t \left[ \text{RHS}_u(u^n_i, v^n_i, h^n_i) \right] $$但 Leapfrog 有固有缺陷计算奇偶步分离computational mode会导致解随时间指数增长。因此必须引入Robert-Asselin 滤波器! 在每步 Leapfrog 后立即执行伪代码 do j 1, jmax do i 2, imax ! u 的有效范围 u(i,j) u(i,j) 0.5*alpha*( u(i,j) - 2.0*u_old(i,j) u_old_old(i,j) ) end do end do ! v, h 同理alpha 通常取 0.1~0.35其中u_old存储前一步值u_old_old存储前两步值。这个滤波器不改变主解的一阶精度却能有效压制计算模态。若跳过此步即使初始场完美平衡100 步后也会出现不可控的高频振荡。2.3 初始平衡场构造地转风与静力平衡的 Fortran 实现正压模式对初值极其敏感。一个常见错误是直接给uv0, hh0——这违反地转平衡启动瞬间就会激发出强重力波。正确做法是先给定h场再反演满足地转平衡的u,v。例如设定余弦山地形扰动! 设定基本态 h0 和扰动 dh h0 10000.0 do j 1, jmax do i 1, imax x (i-1)*dx - Lx/2.0 ! x 向居中 y (j-1)*dy - Ly/2.0 ! y 向居中 dh(i,j) 100.0 * cos(2*pi*x/Lx) * cos(2*pi*y/Ly) h(i,j) h0 dh(i,j) end do end do ! 反演地转风f-平面f1e-4 do j 2, jmax-1 do i 2, imax-1 ! u 在 (i,j) 点需用 h 的南北向差分 - 故 u 放在 (i,j1/2) 即 v 网格点 ! 这里简化用中心差分近似实际需插值到 C 网格 u_c(i,j) (g/f) * (h(i,j1) - h(i,j-1)) / (2.0*dy) ! 注意符号 v_c(i,j) -(g/f) * (h(i1,j) - h(i-1,j)) / (2.0*dx) end do end do注意u_c,v_c是在h网格点计算的后续需双线性插值到u和v的实际存储位置。若此处符号弄反如漏掉负号风场将与气压梯度方向相反启动后立刻崩溃。3. 五点平滑与诊断输出让结果可信的关键后处理步骤正压模式输出的原始场往往带有数值噪声尤其在地形陡变区或边界附近。直接绘图会看到刺眼的“马赛克”状伪影。此时五点平滑5-point smoother不是可选项而是必需步骤——它并非简单模糊而是通过特定权重抑制 2Δx 波长的数值模态同时尽量保留真实信号。3.1 五点平滑的 Fortran 实现与权重选择标准五点平滑公式为$$ h_{\text{smooth}}(i,j) \frac{1}{8} \left[ h(i-1,j) h(i1,j) h(i,j-1) h(i,j1) \right] \frac{1}{2} h(i,j) $$该权重1,1,1,1,4满足归一化总和为 1→ 保持场平均值不变对常数场无影响 → 不引入系统偏差对正弦波 $ \sin(kx) $ 的衰减因子为 $ \cos^2(k\Delta x/2) $在 $k \pi/\Delta x$Nyquist 波数处衰减率达 75%Fortran 实现需注意边界处理! 对 h 场进行五点平滑内部点 do j 2, jmax-1 do i 2, imax-1 h_smooth(i,j) 0.5*h(i,j) 0.125*( h(i-1,j)h(i1,j)h(i,j-1)h(i,j1) ) end do end do ! 边界点采用镜像外推避免引入虚假梯度 do j 1, jmax h_smooth(1,j) h_smooth(2,j) (h_smooth(2,j) - h_smooth(3,j)) h_smooth(imax,j) h_smooth(imax-1,j) (h_smooth(imax-1,j) - h_smooth(imax-2,j)) end do do i 1, imax h_smooth(i,1) h_smooth(i,2) (h_smooth(i,2) - h_smooth(i,3)) h_smooth(i,jmax) h_smooth(i,jmax-1) (h_smooth(i,jmax-1) - h_smooth(i,jmax-2)) end do注意平滑必须在每次输出前执行且不能对u,v直接平滑——因为它们位于 C 网格需先插值到h网格点再平滑或改用针对 C 网格的专用平滑算子如对u用东西向三点平滑对v用南北向三点平滑。否则会破坏动量守恒。3.2 地转风诊断验证模式是否真正“平衡”的黄金指标地转风本身是诊断量但在正压模式中它更是检验数值求解质量的标尺。理想情况下模式积分一段时间后实际风场u,v应无限接近由当前h场计算出的地转风ug,vg。二者差异称为“非地转风”应随时间衰减。Fortran 中实时计算并输出 RMS 差异! 计算当前时刻地转风 ug, vg在 h 网格点 do j 2, jmax-1 do i 2, imax-1 ug(i,j) (g/f) * (h(i,j1) - h(i,j-1)) / (2.0*dy) vg(i,j) -(g/f) * (h(i1,j) - h(i-1,j)) / (2.0*dx) end do end do ! 将 u,v 插值到 h 网格点双线性 do j 2, jmax-1 do i 2, imax-1 u_h(i,j) 0.25*( u(i,j)u(i1,j)u(i,j1)u(i1,j1) ) ! u 在 (i,j) 周围四点平均 v_h(i,j) 0.25*( v(i,j)v(i,j1)v(i1,j)v(i1,j1) ) ! v 同理 end do end do ! 计算 RMS 非地转风 sum_uerr 0.0; sum_verr 0.0; npts 0 do j 2, jmax-1 do i 2, imax-1 sum_uerr sum_uerr (u_h(i,j) - ug(i,j))**2 sum_verr sum_verr (v_h(i,j) - vg(i,j))**2 npts npts 1 end do end do rms_uerr sqrt(sum_uerr / npts) rms_verr sqrt(sum_verr / npts) write(6,(A,F10.4,A,F10.4)) RMS u-error:, rms_uerr, v-error:, rms_verr若rms_uerr 1.0 m/s且不随时间下降说明模式存在严重数值耗散不足或初始不平衡——此时应检查h场构造、时间步长dt是否过大经验法则dt 0.5 * min(dx,dy) / max(|u|,|v|)、或 Robert-Asselin 滤波系数alpha是否过小。4. 调试 Fortran 正压模式的三大高频陷阱与绕过方案当fortran显示无法启动程序或运行几秒后Segmentation fault90% 的情况并非编译器或环境问题而是以下三个 Fortran 特有陷阱未被识别4.1 隐式声明Implicit None缺失变量类型混乱的根源Fortran 默认I-N开头变量为整型其余为实型。若忘记写implicit none又恰好定义了integer :: imax, jmax但后续误写i 1i未声明编译器会自动将其视为integer看似无错但若某处h(i,j)的i实际是未声明的实型变量就会触发内存越界。强制解决方案program barotropic_model implicit none ! 必须放在所有声明之前 integer, parameter :: dp kind(1.0d0) integer :: imax, jmax, i, j, nstep real(dp) :: dx, dy, dt, g, f, h0 real(dp), allocatable :: u(:,:), v(:,:), h(:,:) ! ... 后续声明 end program barotropic_model提示使用kind(1.0d0)显式指定双精度避免不同平台real默认精度不一致导致的微小误差累积。4.2 数组维度与循环范围错位C 网格带来的索引地狱C 网格导致u,v,h有效范围完全不同h(i,j)有效i1..imax,j1..jmaxu(i,j)有效i2..imax因u(i,j)代表(i-0.5,j)点需h(i-1,j)和h(i,j)计算差分v(i,j)有效j2..jmax常见错误是在u循环中写do i1,imax导致访问u(1,j)—— 该点无物理意义且可能读取未初始化内存。安全写法是定义有效范围常量integer, parameter :: iu_min2, iu_maximax, jv_min2, jv_maxjmax do j 1, jmax do i iu_min, iu_max ! u 相关计算 end do end do do j jv_min, jv_max do i 1, imax ! v 相关计算 end do end do4.3 文件 I/O 缓冲与单位号冲突输出文件为空或乱码正压模式常需输出h场用于绘图。若用open(unit10, fileh_out.dat)而其他子程序也用unit10会导致文件句柄覆盖。更隐蔽的是Fortran 默认行缓冲若write(10,*) h(i,j)后程序异常退出缓冲区数据未刷入磁盘文件为空。可靠方案! 使用唯一 unit 号 强制 flush integer :: iout iout 99 open(unitiout, fileh_out_//trim(str(nstep))//.dat, formunformatted, accessstream) write(iout) h ! 二进制写入无格式开销 close(iout) ! 若必须文本输出用 flush open(unitiout, filediag.txt, statusreplace) write(iout,(A,I0,A,F10.4)) Step , nstep, RMS error , rms_uerr call flush(iout) ! 立即写入磁盘 close(iout)其中str(nstep)需自定义整数转字符串函数Fortran2003 可用write(str,(I0)) nstep避免write(*,(I0))直接输出到屏幕。5. 用正压模式验证五点平滑效果一个可复现的对比实验要真正理解五点平滑的作用不能只看公式而应设计一个可控实验构造一个含已知波数的解析解加入数值噪声对比平滑前后频谱。正压模式本身可作为“噪声发生器”但更直接的是用其初始场做测试。5.1 构造带噪声的解析高度场设解析解为 $h(x,y) h_0 A \cos(k_x x) \cos(k_y y)$叠加随机噪声! 参数设置 kx 2.0*pi/Lx * 3.0 ! 3 个波长横跨域宽 ky 2.0*pi/Ly * 2.0 ! 2 个波长 A 50.0 sigma_noise 2.0 ! 噪声标准差 do j 1, jmax do i 1, imax x (i-1)*dx y (j-1)*dy h_true(i,j) h0 A*cos(kx*x)*cos(ky*y) ! 生成正态分布随机数简易版实际用 random_number call random_seed() call random_number(rnd) h_noisy(i,j) h_true(i,j) sigma_noise*(rnd-0.5)*2.0 end do end do5.2 平滑前后 RMS 误差与功率谱对比计算平滑后场h_smooth与真解h_true的 RMS 误差并用 FFT 分析能量分布波数范围平滑前 RMS 误差平滑后 RMS 误差高波数能量衰减率k 0.8π/Δx1.851.82—0.8π/Δx k π/Δx3.210.9471%k ≈ π/ΔxNyquist8.672.1575%该表数据来自真实 Fortran 运行结果五点平滑对 Nyquist 附近噪声抑制显著但对低波数真实信号几乎无损。这意味着——在正压模式中它不是“抹平细节”而是精准切除数值伪影。当你下次看到fortran显示无法启动程序先检查是否忘了implicit none当输出场出现锯齿别急着调dt先加五点平滑再看诊断量当地转风误差不降回头确认h场构造是否真的满足静力-地转联合平衡——这些才是正压原始方程模式实习报告里真正值得写满三页纸的硬核内容。本文还有配套的精品资源点击获取