简介：POM普林斯顿海洋模型是一套面向海洋数值模拟初学者与科研人员的建模工具，可在f77、f90、g77、pgf77、CVF等编译器中直接运行，用于建立河道水流与河口区域的水动力模型。资源包共21个文件，约658KB，包含f、c源码文件、out输出结果、n与dat数据文件、runpom2k运行脚本、pdf中文说明书及zip、tar打包的辅助工具，覆盖从代码编译到结果查看的完整流程。POM2K提供两个实例：河道模型参数已内置于代码中，河口模型则附带原始网格文件，其输入文件由GRID-DATA/GRID.f生成，用户可据此设置自己的网格并插值地形、温盐、表面风场及热能量等数据，生成初始条件与强迫条件。目前已有997人学习下载，适合需要快速上手海洋模型配置、理解网格插值与输入文件生成流程的读者参考。
1. 为什么一个四十年前的海洋模型，今天还在被反复重写
如果你翻过任何一份区域海洋数值模拟的代码库，大概率会撞见三个字母：POM。它诞生于上世纪八十年代，比很多读者的年龄都大，却依然活跃在河口、陆架、近海环流的研究与工程代码里。反直觉的地方在于：它并不是精度最高的模型，也不是功能最全的模型，甚至它的原始代码风格在今天看来相当“上古”——但偏偏是它，成了无数后来者（ROMS、FVCOM 等）在结构设计上的参照系。原因不复杂：POM 用一套极其克制的数学离散方案，把“自由表面、斜压、湍流闭合”这三件海洋模拟里最麻烦的事，塞进了一个可以被人完整读懂的代码框架里。搞懂 POM 的底层机制，等于拿到了一把钥匙——你能看懂为什么后来的模型要那样切网格、那样分步积分、那样处理垂向混合。这份资源正是围绕“建立模型”这个目标，把 POM 的骨架拆开给你看。
2. 控制方程与离散：POM 到底在解什么
2.1 原始方程组：海洋模拟绕不开的六个未知量
任何海洋模型的第一步都是把海水当成连续介质，写出动量、连续、温度、盐度守恒方程。POM 采用的是静力近似与 Boussinesq 近似下的原始方程组，未知量包括水平流速 u、v，垂向流速 w，自由表面高度 η，位温 T 和盐度 S。静力近似意味着垂向动量方程退化为压力梯度与浮力平衡，这一步砍掉了垂向声波与快速重力波，代价是不能模拟对流尺度的垂向加速度，收益是时间步长可以放大几个量级。Boussinesq 近似则把密度变化只保留在浮力项里，其余位置用参考密度替代，避免了解一个变密度不可压流体的复杂状态方程。
常见的实现思路是：先把三维方程沿垂向做 sigma 坐标变换，把随海底起伏的物理域映射成规则的计算域。sigma 坐标的定义是 σ = (z - η)/(H + η)，从海底的 -1 到表面的 0。这样表面和海底都变成坐标面，边界条件处理起来干净很多，代价是压力梯度项在陡坡处会出现“sigma 坐标压力梯度误差”，这是 POM 类模型最著名的坑之一，后面避坑章节会专门讲。
2.2 模态分裂：为什么要把正压和斜压拆开算
POM 最核心的设计决策是模态分裂（mode splitting）。自由表面重力波传播速度是 √(gH)，远快于内波和平均流的平流速度。如果用一个统一的时间步长同时算两者，要么步长被正压波压得极小，要么斜压过程精度不够。POM 的做法是：把解拆成“正压模态”（垂向平均的流速与水位）和“斜压模态”（偏离垂向平均的部分），正压用短步长 Δt_p 积分，斜压用长步长 Δt_b 积分，两者步长比通常在 10 到 50 之间。
c POM 主循环的模态分裂骨架（示意） do 100 n=1,nstep c 正压子循环：用短步长推进水位和垂向平均流速 do 200 nps=1,np call advave ! 正压平流项 call extel ! 外模式水位与输运 200 continue c 斜压步：用长步长推进三维温度、盐度、湍流闭合 call profq ! 垂向扩散系数 call advq ! 斜压平流 call dens ! 状态方程更新密度 call baropg ! 斜压压力梯度 100 continue这段骨架里，extel负责外模式（正压）的水位与输运更新，baropg负责把斜压密度场产生的压力梯度反馈给三维动量方程。两个模态之间通过垂向平均流速和斜压压力梯度双向耦合。参数上，np是每个斜压步里正压子循环的次数，它由 CFL 条件决定：正压波速乘以 Δt_p 必须小于网格间距。如果np设得太小，正压模态会数值不稳定；设得太大，计算时间线性增长。我一般会这么理解：模态分裂本质上是把“快过程”和“慢过程”在时间尺度上解耦，用不同的积分精度去匹配各自的物理特征时间，这是 POM 能在当年有限算力下跑出可用结果的关键。
2.3 垂向湍流闭合：Mellor-Yamada 2.5 阶方案在算什么
POM 默认搭载的是 Mellor-Yamada 2.5 阶湍流闭合模型。它不直接解湍流脉动，而是解湍动能 q² 和湍流宏观尺度 l 的输运方程，再用它们构造垂向涡黏系数 K_M 和涡扩散系数 K_H。2.5 阶的意思是：保留湍动能和尺度两个预报量，但对方程中的三阶矩项做参数化截断。常见的实现思路是，q² 方程包含产生项、耗散项、扩散项和浮力项，l 方程则通过一个经验性的壁面函数在边界附近约束尺度。
c 计算垂向涡黏系数（示意） do k=1,kbm1 q2 = max(q2b(k), q2min) ! 限制最小湍动能，防止除零 l = max(lb(k), lmin) ! 限制最小尺度 km(k) = l * sqrt(q2) * sm ! sm 为稳定度函数 kh(k) = l * sqrt(q2) * sh end doq2min和lmin是数值保护下限，防止在层结极强或流速极弱时湍动能趋近于零导致系数矩阵奇异。sm和sh是 Richardson 数的函数，层结越稳定，它们越小，垂向混合越弱。这套闭合方案的好处是物理可解释性强，参数少；坏处是对流边界层和强风搅拌的刻画偏弱，这也是后来很多模型改用 k-ε 或 k-ω 方案的原因。但在 POM 的框架里，Mellor-Yamada 与 sigma 坐标、模态分裂是配套设计的，换闭合方案往往需要重新调步长和滤波参数。
3. 网格、边界与数值格式：模型能跑起来的前提
3.1 Arakawa C 网格：为什么速度不放在格点上
POM 水平方向采用 Arakawa C 网格：水位、温度、盐度放在网格中心，u 速度放在东西向格点之间，v 速度放在南北向格点之间。这种交错布置的好处是，压力梯度项和散度项可以用相邻格点的差分直接写出，不需要插值，且能有效抑制棋盘式压力振荡。代价是动量平流项和科氏力项需要做四点平均，代码里会出现大量0.25*(u(i,j)+u(i+1,j)+...)形式的平均操作。
垂向采用 sigma 坐标后，每一层的水深随位置变化，所以垂向差分也要乘上对应的层厚因子。POM 里用dz数组存储各层厚度，动量方程中的垂向扩散项写成(km(k+1)*(u(k+1)-u(k))/dz(k+1) - km(k)*(u(k)-u(k-1))/dz(k))/dz(k)这类形式。参数上，kbm1是垂向层数减一，层数越多，对温跃层的刻画越细，但计算量和垂向 CFL 约束也越紧。
3.2 边界条件：开边界、海底应力和表面风应力
开边界是区域海洋模型最容易翻车的地方。POM 支持几种开边界处理：辐射边界条件、指定水位边界、以及嵌套边界。辐射条件的形式是 ∂φ/∂t + c ∂φ/∂n = 0，其中 c 是外传波速，通常取正压波速或流速。如果 c 取错符号，边界会反射而不是吸收，导致能量在区域内堆积，模拟结果出现虚假振荡。
海底边界条件通过底摩擦应力与流速的平方关系给出，拖曳系数 C_d 通常取 0.0025 左右，但在浅水区需要根据底质类型调整。表面边界条件则接收风应力、热通量和淡水通量。风应力的方向决定上层埃克曼输运的方向，热通量决定层结强度，淡水通量影响盐度与密度。常见的实现思路是：把这些通量作为强迫场按时间插值输入，时间分辨率至少要到日平均，否则会漏掉天气尺度的风事件。
3.3 时间积分格式： leapfrog 与 Asselin 滤波
POM 的时间积分主格式是 leapfrog（蛙跳）格式，它对线性平流有很好的守恒性，但存在计算模态（computational mode）问题：奇数步和偶数步的解会逐渐分离，产生锯齿状振荡。POM 用 Asselin 时间滤波来抑制这个模态，滤波系数通常取 0.05 到 0.1。
c Asselin 滤波（示意） do k=1,kb uf(k) = u(k) + 0.5*alpha*(uold(k) - 2.0*u(k) + unew(k)) end doalpha是滤波系数，越大抑制越强，但也会削弱物理信号的振幅。这个参数没有普适最优值，需要根据模拟区域和关注的时间尺度试。如果模拟的是潮汐，alpha取小一点保留潮汐振幅；如果模拟的是长期环流，可以取大一点换稳定性。
4. 避坑与排查：POM 建模中最容易翻车的五件事
4.1 现象：模拟跑了几十步就发散，流速指数增长
原因通常是正压 CFL 条件被违反。正压波速 √(gH) 在深水区很大，如果 Δt_p 没有相应缩小，或者np与 Δt_b 的比值设置不当，外模式会先崩。解决方法是先算一遍各网格点的 √(gH)·Δt_p/Δx，确保最大值小于 1，然后反推np。另外检查开边界处是否有异常大的水深输入。
4.2 现象：温度场出现棋盘状振荡
原因多半是 Arakawa C 网格上的平流格式选择不当，或者时间滤波系数太小。POM 提供了中心差分和迎风差分两种平流选项，中心差分精度高但在强梯度处容易振荡，迎风差分稳定但数值扩散大。解决方法是先切到迎风差分确认振荡消失，再逐步调回中心差分并增大 Asselin 系数。
4.3 现象：陡坡区域出现虚假的强垂向混合
这是 sigma 坐标压力梯度误差的典型表现。在海底坡度大的地方，压力梯度项的两个大项相减，舍入误差被放大，产生虚假的密度梯度，进而驱动虚假的垂向混合。解决方法是采用地形跟随坐标的压力梯度修正算法，或者对地形做平滑，把相邻网格的水深比控制在 1.2 以内。
4.4 现象：开边界附近盐度或温度异常堆积
原因是辐射边界条件的波速符号或大小不对，导致边界反射。解决方法是检查边界处的外传波速是否与内部流速方向一致，必要时改用指定边界值或嵌套方式。如果模拟区域有河流输入，还要确认淡水通量的位置和量级是否与边界条件协调。
4.5 现象：湍流闭合算出的 K_M 在表层异常大
原因通常是风应力强迫过强，或者 q² 方程的边界条件设置不当。Mellor-Yamada 方案在表层需要给定 q² 和 l 的边界值，如果 l 的壁面函数参数取错，表层混合层会被人为加深。解决方法是检查q2min、lmin和壁面函数的系数，并与实测的混合层深度做对比。
5. 从 POM 到现代模型：设计取舍与验证思路
POM 的设计取舍在今天看来非常清晰：它用 sigma 坐标换取了海底边界的贴合，用模态分裂换取了正压波的积分效率，用 Mellor-Yamada 闭合换取了物理可解释性。这三者互相咬合，构成了一个自洽的框架。但代价也很明显：sigma 坐标在陡坡处的压力梯度误差、模态分裂对步长比的敏感性、以及 2.5 阶闭合对强对流和波浪破碎的无力。后来的模型有的改用 z 坐标或混合坐标来规避压力梯度误差，有的用更高级的湍流闭合，有的用非结构网格来适应复杂岸线——但如果你仔细看它们的代码结构，模态分裂的思想、交错网格的布置、时间滤波的策略，很多都能追溯到 POM。
验证一个 POM 模型是否可信，我一般会走三步。第一步是保守示踪剂测试：给一个初始浓度场，不加任何强迫，看它是否在预期时间内均匀混合，用来检查数值扩散和边界处理。第二步是潮汐验证：用实测潮位数据对比模拟水位，看振幅和相位误差是否在可接受范围，这一步能暴露正压模态的问题。第三步是温盐剖面验证：用 CTD 或浮标数据对比垂向结构，看混合层深度和温跃层位置是否合理，这一步能暴露湍流闭合和斜压压力梯度的问题。
| 验证项 | 关注指标 | 常见误差来源 |
|---|---|---|
| 保守示踪剂 | 浓度方差衰减速率 | 数值扩散、时间滤波系数 |
| 潮汐水位 | 振幅误差、相位差 | 正压步长、开边界条件 |
| 温盐剖面 | 混合层深度、温跃层梯度 | 湍流闭合、sigma 压力梯度 |
| 流速剖面 | 垂向结构、底摩擦 | 底拖曳系数、垂向分辨率 |
从那以后我每次拿到一个新的 POM 配置，都会先跑一遍保守示踪剂测试，再对照至少一个潮位站的数据，最后才去看温盐结果。这个习惯帮我省掉了大量“看起来像物理、其实是数值”的误判。希望帮到你。
本文还有配套的精品资源，点击获取
