企业官网建设流程全解析

简介：POM普林斯顿海洋模型是一套面向海洋数值模拟初学者与科研人员的开源建模工具，支持在f77、f90、g77、pgf77、CVF等多种Fortran编译器中直接运行，可用于河道水流与河口等典型场景的建模实践。资源包共21个文件，约658KB，包含f、c源码文件、out输出结果、n与dat数据文件、runpom2k运行脚本，以及POM中文说明书PDF和pomviz可视化压缩包，覆盖从编译、运行到结果查看的完整流程。POM2K提供两个可直接上手的实例：河道模型参数已内置于代码中，河口模型则附带原始网格文件，其输入文件由GRID-DATA/GRID.f生成，用户可据此设置自己的网格并插值地形、温盐、表面风场与热通量等资料，生成初始条件与强迫条件。目前已有997人学习下载，适合希望快速理解POM建模流程、掌握网格与数据插值方法的读者参考。

1. 普林斯顿海洋模型：为什么一个四十年前的代码，今天还在跑

如果你翻开源海洋模式的谱系，会发现一个反直觉的事实：今天很多业务化预报系统、区域海洋再分析产品、甚至耦合气候模式里的海洋分量，往上追溯都能碰到同一份代码——普林斯顿海洋模型（POM）。它诞生于上世纪八十年代末，用 Fortran 77 写成，没有花哨的并行框架，没有自动微分，连网格都是最简单的 Arakawa C 交错网格。可就是这样一份"老古董"，至今仍被大量研究组当作区域海洋模拟的起点。

这件事值得琢磨。海洋模式不是越新越好吗？为什么一个没有 GPU 加速、没有非结构网格、连湍流闭合都停留在 2.5 阶的老模型，还能在近岸、河口、陆架这些最复杂的海域里稳定工作？答案藏在它的设计哲学里：POM 用一套"自由表面 + 模态分裂 + 二阶湍流闭合"的组合，把三维海洋的斜压演化问题拆成了两个时间尺度上可解的子问题。搞懂这套拆分逻辑，你就能理解为什么它能在几十年的迭代里活下来，也能理解今天那些更复杂的模式在继承什么、又在改什么。

这份资源是一套完整的 POM 建模资料，从控制方程到数值离散，从编译配置到算例搭建都有覆盖。下面不按"怎么装怎么跑"来讲，而是把它的底层机制一层层拆开——先看它怎么把物理问题变成可计算的数学结构，再看这些结构在源码里长什么样，最后聊它在什么场景下会翻车、以及为什么这些取舍至今仍有参考价值。

2. 模态分裂：POM 为什么要把海洋拆成快慢两层

2.1 正压模态与斜压模态的时间尺度差

海洋里有两种运动，速度差了两个数量级。一种是表面重力波和正压潮汐，传播速度是 √(gH)，取水深 100 米，大约是 31 m/s；另一种是内波和斜压环流，速度由约化重力 g' = g·Δρ/ρ 决定，典型值只有 1～3 m/s。如果用一个显式时间积分同时算这两种运动，时间步长会被最快的表面波卡死——CFL 条件要求 Δt < Δx/√(gH)，取 Δx = 1 km，Δt 只有几十秒。可你真正关心的斜压过程，时间尺度是小时到天。用几十秒的步长去算几天的环流，计算量直接爆炸。

POM 的解法是模态分裂（mode splitting）：把三维速度分解成"垂直平均 + 偏差"两部分。垂直平均后的正压模态用短步长积分，只算表面高度和正压输运；斜压模态用长步长积分，算三维温度、盐度和偏差速度。两个模态通过表面高度和正压输运耦合。这样短步长只花在二维问题上，长步长承担三维计算，总计算量降下来一个量级。

常见的实现思路是：正压模态用显式格式，斜压模态里对平流项用显式、对垂直扩散项用隐式，避免垂直方向的扩散限制步长。这个"显隐混合"的选择，是 POM 能在近岸高分辨率下稳定运行的关键。

2.2 分裂步长的约束与耦合方式

模态分裂不是随便拆的，它有一个硬约束：正压模态的积分必须"追上"斜压模态的时间。具体说，如果斜压步长是 Δt，正压步长是 Δt/n，那么一个斜压步里要跑 n 次正压积分。n 的取值由 CFL 条件决定：

n ≥ ceil( Δt · √(gH_max) / Δx_min )

这个 n 不能太小，否则正压模态不稳定；也不能太大，否则正压部分的计算量反超斜压部分，分裂就失去意义。实际配置里，n 通常在 10～50 之间，取决于网格分辨率和最大水深。

耦合方式上，POM 用的是"预测-校正"思路：斜压步开始时，用当前表面高度算斜压速度的垂直积分，得到正压输运的初始场；正压步在 n 个小步里推进表面高度和输运；斜压步结束时，用新的表面高度更新三维速度。这个过程中，表面高度是唯一的耦合变量，正压模态不碰温度盐度，斜压模态不碰快波——职责边界非常清晰。

注意：模态分裂的稳定性依赖于正压模态的 CFL 满足，如果网格里有特别浅的点（比如滩涂），√(gH) 变小，n 反而可以取小；真正卡步长的是最深的地方。配置时按最大水深算 n，不要按平均水深。

2.3 源码里的模态循环结构

在 POM 的主程序里，模态分裂体现为两层循环。下面是一段简化后的结构，保留了关键调用关系：

! 斜压步循环：advance 三维温盐和偏差速度 do k = 1, nbc ! nbc 是斜压步数 ! 计算斜压项：平流、水平扩散、垂直扩散 call advct ! 温盐平流 call profq ! 垂直扩散系数 call advq ! 温盐垂直扩散 call baropg ! 斜压压强梯度力 ! 正压子循环：n 个小步推进表面高度和输运 do i = 1, n call bcalc ! 算正压速度 call pgrad ! 正压压强梯度 call advave ! 正压平流 call extel ! 推进表面高度 end do ! 用新的表面高度更新三维速度 call vertvl ! 垂直速度 call dens ! 密度 end do

这段代码里，bcalc和extel是正压模态的核心，advct、profq、advq是斜压模态的核心。baropg算的是斜压压强梯度力，它把密度场的信息传给动量方程。整个结构的精妙之处在于：正压子循环里没有任何三维数组的垂直运算，全是二维场操作，缓存友好，向量化容易。这也是 POM 在当年硬件上能跑得动的原因。

3. 湍流闭合：2.5 阶 Mellor-Yamada 到底在算什么

3.1 为什么海洋模式必须做湍流闭合

海洋模式的控制方程是雷诺平均后的 Navier-Stokes 方程，平均过程会多出一个雷诺应力项，代表湍流对平均流的输运。这一项不闭合，方程就没法解。所谓"湍流闭合"，就是用一个参数化的方式，把雷诺应力和大尺度变量（速度剪切、温度梯度、浮力频率）联系起来。

最简单的做法是常系数扩散，但海洋里的垂直混合系数变化范围极大：表层风混合可以到 10⁻² m²/s，温跃层里只有 10⁻⁵ m²/s，差三个量级。用常数等于放弃治疗。POM 选的是 Mellor-Yamada 2.5 阶闭合，它不解完整的湍流应力方程，而是解湍动能（TKE）的输运方程，再用 TKE 和长度尺度算混合系数。这个"2.5 阶"的意思是：保留了湍动能的预报方程（一阶矩），但应力和通量用代数关系闭合，介于 1 阶和 2 阶之间。

3.2 TKE 方程与稳定函数

Mellor-Yamada 2.5 阶的核心是湍动能 q²/2 的预报方程：

∂(q²/2)/∂t + u·∇(q²/2) = ∂/∂z[ Kq · ∂(q²/2)/∂z ] + P + B - ε

其中 P 是剪切产生项，B 是浮力产生/耗散项，ε 是耗散率。POM 里把 q² 记作q2，长度尺度记作q2l，两个量都做预报。混合系数由稳定函数算：

Km = l · q · Sm ! 动量混合系数 Kh = l · q · Sh ! 温盐混合系数

Sm 和 Sh 是稳定函数，依赖于理查森数 Ri = N²/(∂u/∂z)²。当 Ri 大（层结强），Sm、Sh 变小，混合被抑制；当 Ri 小（剪切强），混合增强。这个反馈机制让模式能自动区分温跃层和混合层，不需要人为设阈值。

源码里，profq子程序负责算 q2 和 q2l 的垂直扩散系数，advq负责推进这两个量。稳定函数的计算在prodq和dissq里，涉及一堆经验常数（比如 B1 = 16.6，E1 = 1.8，E2 = 1.33）。这些常数来自实验室和观测拟合，不是从第一性原理推出来的——这也是 2.5 阶闭合被诟病的地方：它在大洋深处表现不错，但在强层结的河口锋面附近经常偏弱。

3.3 源码里的闭合实现与参数含义

下面这段代码展示了 q2 和 q2l 的推进逻辑，保留了关键计算：

! 计算湍动能 q2 和 q2l 的垂直扩散 do k = 2, kb-1 ! 稳定函数：依赖理查森数 r = gh(k) * (b(k) - b(k+1)) / ((u(k)-u(k+1))**2 + (v(k)-v(k+1))**2 + eps) ! Sm, Sh 由 r 查表或公式算 sm = ...; sh = ... ! 混合系数 km(k) = l(k) * sqrt(q2(k)) * sm kh(k) = l(k) * sqrt(q2(k)) * sh ! 扩散系数 kq(k) = l(k) * sqrt(q2(k)) * 0.2 ! 0.2 是经验常数 end do ! 推进 q2 和 q2l do k = 2, kb-1 q2(k) = q2(k) + dt * (diff(k) + prod(k) - diss(k)) q2l(k) = q2l(k) + dt * (diff_l(k) + prod_l(k) - diss_l(k)) end do

gh是重力加速度乘层厚，b是浮力，r就是局部理查森数。km、kh、kq分别是动量、温盐、湍动能的垂直混合系数。注意kq里的 0.2 是经验常数，对应 Mellor-Yamada 原始文献里的 S_q。这些常数在 POM 的不同版本里有微调，但量级不变。

提示：如果你在河口或强潮海域跑 POM，发现混合层太浅或温跃层太锐，先别急着改网格，去检查profq里的稳定函数。很多"模式跑不通"的问题，根源是 Ri 在强剪切区算出了负值，导致 Sm 爆炸。常见的做法是给 Ri 加一个下限，或者在计算剪切时加一个速度下限。

4. 数值离散：Arakawa C 网格与自由表面怎么配合

4.1 C 网格的交错布局

POM 用的是 Arakawa C 网格：速度分量放在网格单元的面心上，标量（温度、盐度、表面高度）放在单元中心。这种交错方式的好处是，压强梯度力和散度算起来不会出现棋盘格振荡——如果速度和标量放在同一位置，中心差分算出来的压强梯度在网格尺度上会解耦，出现两个相邻格点压强交替高低的假解。C 网格把速度错开半个格点，差分算子天然耦合，避开了这个问题。

具体布局是：u 速度在东西面，v 速度在南北面，w 在垂直层界面，温度盐度在单元中心。表面高度也在中心。这种布局下，连续方程的离散是精确的——每个单元的体积通量用面心速度算，进出平衡直接对应中心的高度变化。

4.2 自由表面与刚盖近似的取舍

早期海洋模式多用刚盖近似：表面高度不预报，而是解一个椭圆方程让垂直积分散度为零。这样做的好处是滤掉了表面重力波，时间步长可以放大；坏处是要解 Poisson 方程，并行效率低，而且不能模拟潮汐。

POM 选了自由表面：表面高度做预报，允许表面起伏。代价是必须处理快波，所以才有模态分裂。这个选择在当年是有争议的——自由表面的计算量更大，但 POM 的作者赌的是"潮汐和风暴潮是近岸海洋的核心过程，不能滤掉"。后来的发展证明这个赌注是对的：区域海洋模式如果不算潮汐，近岸的混合、输运、锋面位置全都不对。

源码里，表面高度的推进在extel子程序：

! 推进表面高度：基于正压输运的散度 do j = 2, jm-1 do i = 2, im-1 ! 东西向输运差 + 南北向输运差 div = (d(i+1,j)*uab(i+1,j) - d(i,j)*uab(i,j)) / dx + (d(i,j+1)*vab(i,j+1) - d(i,j)*vab(i,j)) / dy ! 表面高度更新 el(i,j) = el(i,j) - dt * div end do end do

d是水深，uab、vab是垂直积分的正压输运，el是表面高度。这段代码就是连续方程的离散：高度变化等于输运的净流入。注意d*uab是单宽输运，除以dx得到散度。这个格式是显式的，所以正压步长受 CFL 限制。

4.3 垂直坐标与地形跟随

POM 用 sigma 坐标：垂直方向归一化到 [0,1]，0 是表面，1 是底。这样自由表面和海底地形都能被精确贴合，近岸的底边界层不需要阶梯近似。代价是 sigma 坐标在陡坡上会产生压强梯度力误差——两个相邻格点的密度差乘以坡度，可能比真实的水平压强梯度还大，导致虚假流动。

常见的缓解办法是算压强梯度力时减去一个参考密度剖面，只算偏差部分。POM 的baropg子程序里就有这个处理：先用平均密度算一个背景压强，再从总压强里减掉。这个技巧在陡坡海域是必须的，否则模式会自己"造"出沿岸流。

5. 避坑与排查：POM 建模里最容易翻车的几个地方

5.1 模式跑飞：先查 CFL 还是先查闭合

现象：模式跑了几百步后速度突然发散，表面高度出现量级异常。

原因：九成是 CFL 条件被破坏。要么是正压子步数 n 设小了，要么是网格里有特别浅的点导致 √(gH) 算错，要么是初始场里有速度突变。

解决：先打印最大正压波速和实际步长，算一下 CFL 数。如果大于 1，把 n 调大。如果 n 已经很大还是发散，检查地形数据里有没有负水深或零水深——POM 对零水深没有保护，除零会直接产生 NaN。

5.2 温盐场出现负值：平流格式的锅

现象：跑一段时间后，某个角落的盐度变成负数，或者温度超过物理上限。

原因：POM 默认用的是中心差分平流，在强梯度区会产生过冲。河口锋面、温跃层这些地方梯度大，最容易出问题。

解决：换成迎风或 TVD 格式。POM 的advct里有选项，但默认不一定开。如果不想改格式，至少加一个限幅器，把温盐限制在物理范围内。这个坑在近岸模拟里非常常见，尤其是淡水输入强的河口。

5.3 表面高度漂移：正压模态没收敛

现象：模式跑久了，区域平均表面高度缓慢上升或下降，没有回到零。

原因：正压模态的输运和斜压模态的垂直积分不一致，导致质量不守恒。常见于开边界处理不当，或者正压子步数 n 不是整数倍关系。

解决：检查开边界的流量是否闭合，表面高度边界条件是否和输运匹配。另外，n 必须是整数，且斜压步和正压步的时间要对齐。如果用了自适应步长，确保正压步长始终满足 CFL。

5.4 编译通过但结果全零：初始化顺序错了

现象：模式正常退出，但所有输出变量都是零或常数。

原因：POM 的初始化有严格顺序：先读地形，再算 sigma 坐标的层厚，再初始化温盐场，最后算密度和压强。如果顺序错了，比如在算层厚之前就调用了密度计算，得到的全是零。

解决：对照主程序的调用顺序，确认init阶段的子程序依赖。特别是dens和baropg，它们依赖d、zz、dd这些几何量，必须在几何初始化之后调用。

5.5 并行效率低：模态分裂的通信瓶颈

现象：用 MPI 并行后，加速比远低于预期，甚至比单核还慢。

原因：正压子循环里每一步都要做全局的输运交换和表面高度更新，通信频率是斜压步的 n 倍。如果区域划分不合理，通信量会吃掉所有计算收益。

解决：把正压子循环里的通信合并，减少同步点。另外，区域划分尽量让每个进程的网格块接近正方形，减少边界面积。如果 n 很大，考虑把正压模态也做区域分解，而不是全局求解。

6. 从 POM 的取舍看区域海洋建模的底层逻辑

POM 的设计里有一组贯穿始终的取舍：用模态分裂换计算效率，用 sigma 坐标换地形贴合，用 2.5 阶闭合换实现简单。这些选择在当年是工程约束下的最优解，今天看依然有参考价值——因为区域海洋建模的核心矛盾没变：你要在有限算力下，同时抓住快过程（潮汐、风暴潮）和慢过程（环流、水团演变），还要让模式在复杂地形上稳定。

模态分裂的代价是耦合误差。正压模态和斜压模态的时间尺度差越大，分裂带来的误差越明显。在强潮海域，表面高度的快速变化会通过压强梯度力反馈给斜压速度，如果正压子步数不够，这个反馈会滞后，导致潮汐相位偏差。我一般会这么理解：模态分裂本质上是一种时间上的多尺度方法，它的精度取决于两个尺度之间的分离程度。分离得好，误差小；分离得不好，就得加子步或者换耦合方式。

sigma 坐标的代价是压强梯度力误差。在坡度超过 0.1 的地方，这个误差可以大到和真实信号同量级。POM 用参考密度剖面减背景的办法缓解，但没有根治。后来的模式有的改用 s-z 混合坐标，有的用地形跟随但加压强梯度力修正，都是在补这个坑。如果你要在陆坡或海山附近做模拟，这一条必须提前评估。

2.5 阶闭合的代价是经验常数的不确定性。B1、E1、E2 这些常数来自特定实验，换一个海域未必适用。而且 2.5 阶闭合假设湍流是局地平衡的，在强非局地输运的场景（比如对流穿透混合层）会失效。常见的改进是加非局地通量项，或者换 k-ε 闭合，但计算量都会上去。

验证一套 POM 配置是否靠谱，我习惯走三步：先跑一个理想算例（比如均匀风应力下的 Ekman 层），看混合层深度和理论值差多少；再跑一个真实海域的短期模拟，和观测的温度盐度剖面对比；最后跑一个长周期算例，检查质量守恒和能量收支。这三步走完，基本能判断配置有没有硬伤。

从那以后我每次搭区域海洋模式，都强制先跑理想算例，不管多急都不跳过。因为模态分裂、sigma 坐标、湍流闭合这三个东西的坑，在理想算例里会暴露得最干净，到了真实海域反而被复杂的边界条件掩盖。这份 POM 资料里把控制方程、离散格式、源码结构都摊开了，适合拿来对着源码逐层验证上面这些机制。希望帮到你。

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