简介POM普林斯顿海洋模型POM2K源码包面向海洋与河口数值模拟方向的初学者及研究者用于快速搭建和测试模型。压缩包共21个文件核心包括可直接在f77、f90、g77、pgf77及CVF等编译器下运行的Fortran源码.f/.c、运行脚本runpom2k、网格数据IC.dat、插值工具GRID-DATA以及POM中文说明书整体仅658KB便于下载后快速部署。模型内置两个可选实例一个河道水流模型参数已在代码中预设可直接运行观察结果另一个河口模型网格文件已直接给出输入文件可由GRID-DATA/GRID.f生成方便对照学习建模流程。已有998人学习下载。借助GRID-DATA插值工具读者可将地形、温盐、表面风场及热通量等实测资料插值到自设网格生成初始条件与强迫条件从而从零建立自己的区域模型并迁移到感兴趣的海域或河口场景。1. 一个海洋模型凭什么让海岸工程沿用三十年接到一个海岸环评项目要在数据稀缺的近岸海域回答“这片湾的潮位和流场到底怎么变”直接搬全球大洋模型分辨率不够自己写水动力方程又明显不现实。这时候 POM 普林斯顿海洋模型就是一个绕不开的选项它用 σ 坐标贴合地形用内外模分裂控制计算效率在近岸、陆架、河口这类地形陡变区域有几十年工程验证。对做水动力、温盐扩散、污染物输运或生态耦合建模的人来说POM 能解决的是“用可接受的算力在海湾尺度跑出物理上可信的三维流场”。本文不是模型手册按我自己从零建立一个 POM 模拟项目的顺序把物理框架、网格生成、初始与边界条件、调参和踩坑串成一条可复现的路径。2. POM 的物理框架与代码骨架读代码之前先搞懂 σ 坐标和内外模2.1 σ 坐标与内外模分裂POM 的两根支柱POM 全称 Princeton Ocean Model最核心的设计是 σ 坐标。常见的 z 坐标模型在深水区层高高、浅水区层高低近岸地形稍微复杂就会在底部分层不完整σ 坐标把垂向坐标固定在表面到海底之间任意水深都保证相同层数。映射关系是 σ (z − η) / (H η)其中 η 是自由面高度H 是静水深。这意味着当水深从 10 米变到 100 米每一层只是被拉伸不会被地形“削掉”。这给近岸建模带来了很大的方便但也埋了一个坑水深变化剧烈的网格点上压力梯度计算容易产生误差后面避坑章节会专门讲。第二个支柱是内外模分裂。POM 把三维水动力方程拆成外模和内模外模是二维的求解自由面和垂向平均流速度快受浅水波 CFL 条件约束时间步长必须小内模是三维的求解温度、盐度和垂向流速结构速度慢时间步长可以放宽。常见做法是内模步长取外模步长的 10 倍左右外模跑 10 步内模跑 1 步通过耦合完成一个完整推进。不理解这个机制直接去调时间步长是许多新手翻车的第一现场——把内模步长当成了唯一旋钮结果外模早就超过了稳定性上限。POM 的垂向湍流闭合用的是 Mellor-Yamada 2.5 阶模型水平扩散采用 Smagorinsky 型参数化不额外引入需要手工调平的水平涡粘系数曲线。这意味着多数情况下默认湍流方案是能用的真正要调的是底摩擦系数、水平分辨率、时间步长和强迫加载方式。2.2 代码包里你只需要盯住这几个文件拿到 POM 源码包之后不要急着逐行读方程。先跑通一个自带案例再替换自己的网格和强迫这是最稳的路径。一个典型的 POM 工程文件构成如下主程序与核心子程序常见命名是 pom 开头的 Fortran 源文件负责时间推进和模式控制。参数输入文件常见命名如 pom.in里面按行读入内外模步长、运行天数、输出间隔、底摩擦参数、σ 分层参数。网格文件通常由专门的网格生成程序产生包含经度、纬度、水深等字段。初始条件和边界条件文件温盐初始场、开边界水位流速强迫。输出文件经典版本多为二进制需要后处理脚本转成可画图的格式。编译和跑通自带案例一般是这样# 1. 检查编译器与依赖库 gfortran --version # 常见依赖为 netCDF Fortran 库经典二进制版则只依赖 Fortran 编译器 # 2. 进入自带案例目录执行编译 make clean make # 3. 跑自带算例看能否正常推进并输出水位场 ./pom pom.in run.log 21 tail -20 run.log逻辑说明先确认编译器存在避免 make 到一半才发现缺依赖make clean 是必要的因为模型源码里常有上一轮编译的中间文件。运行命令把标准输出和错误都重定向到 run.log目的是排查崩溃时能看时间步推进到哪一步。如果 run.log 里出现了 NaN、发散或者水位数值异常先检查编译选项和数据文件路径再看物理参数。参数说明pom.in 是最先要读完整的一个文件常见做法是把每一行参数的默认值先抄一份注释版再逐个改。重点关注运行时长、输出间隔、内模步长、外模步长与内外模步长比。很多 POM 版本的崩溃并非物理参数不合理而是内外模步长比和 CFL 条件不匹配这一点在避坑章节展开。提示不要用 MPI 并行版本去调试物理参数先单进程跑通小网格效率问题永远放在正确性问题之后解决。3. 用深度数据生成 POM 网格正交曲线网格与 σ 分层的落地路径3.1 准备水深数据从地形数据集到网格节点POM 采用正交曲线坐标意味着模型域必须落在“经纬线处处正交”的曲线网格上。这不是随便用矩形网格插值就能交差的长直海岸可以近似用等间距经纬网格但海湾、海峡、岛屿周围需要生成贴合岸线的正交曲线网格。最稳妥的做法是先确定模型域范围和水平分辨率然后把公开地形数据集的水深插值到网格节点上。常见做法是确定模拟范围例如一个海湾东西 40 km、南北 60 km。确定水平网格距近岸工程常见 100 m 到 1 km网格越小越能解析航道和岸线但外模步长随之变小算力成本非线性上升。从地形数据集提取该范围的水深统一到同一个坐标系。插值到模型网格节点输出 POM 网格文件。下面是一段水深插值与网格文件生成的示意脚本import numpy as np import netCDF4 as nc # 读取地形数据集经纬度网格 ds nc.Dataset(bathymetry.nc) lon_src ds.variables[lon][:] lat_src ds.variables[lat][:] dep_src ds.variables[elevation][:] # 注意负值为水下 # 定义模型网格正交曲线网格的 lon/lat 节点 # 这里简化为规则经纬网格真实案例应使用正交网格生成程序的输出 lon_model, lat_model np.meshgrid( np.linspace(121.0, 121.4, 80), np.linspace(30.0, 30.6, 120) ) # 用双线性插值把水深映射到模型网格 from scipy.interpolate import RegularGridInterpolator interp RegularGridInterpolator( (lat_src, lon_src), dep_src, bounds_errorFalse, fill_valueNone ) dep_model interp(np.stack([lat_model.ravel(), lon_model.ravel()], axis-1)) dep_model dep_model.reshape(lon_model.shape) # 转为水下为正的深度 dep_model np.abs(dep_model) # 直接写文本格式的网格文件经度、纬度、水深 grid_out np.stack([lon_model.ravel(), lat_model.ravel(), dep_model.ravel()], axis-1) np.savetxt(pom_grid.dat, grid_out, fmt%.4f %.4f %.1f)逻辑说明地形数据集的海拔负值表示水下插值后取绝对值转成水深。POM 网格文件的具体字段顺序不同版本有差异常见为经度、纬度、水深但有些版本还要求旋转角等参数所以这个输出只能作为中间产物最终字段顺序务必对照代码包里读网格的子程序核对。参数说明插值的网格分辨率决定了地形平滑程度。如果原始地形在 500 m 网格上有 30 m 的突变插值到 1 km 网格后仍然会很陡后续必须做水深平滑否则压力梯度误差会直接表现为近底层的虚假流速。这里有一个工程经验值相邻网格水深变化比一般控制在 0.2 以内超过就要平滑。3.2 正交性检查与水深平滑建模第一个大坑很多人在网格生成之后直接开始跑模型结果在湾口或岛屿附近出现锯齿状的流场。这不是物理问题而是网格正交性不达标。POM 的差分格式假设网格线处处正交实际生成网格很难做到所以每一步都要量化检查。常见检查是逐网格计算两个坐标方向的夹角偏离 90 度的程度用正交角偏差表示工程上常见要求偏差小于 5 度敏感区域小于 2 度。检查脚本如下# 读取网格文件的经度纬度列 lon lon_model.ravel() lat lat_model.ravel() nx, ny 80, 120 # 对每个内部网格点用相邻点差分算出两个切向量再求夹角 ortho_err np.zeros((ny-2, nx-2)) for j in range(1, ny-1): for i in range(1, nx-1): # 沿第一坐标方向的切向量 v1 np.array([lon[j, i1]-lon[j, i-1], lat[j, i1]-lat[j, i-1]]) # 沿第二坐标方向的切向量 v2 np.array([lon[j1, i]-lon[j-1, i], lat[j1, i]-lat[j-1, i]]) # 归一化后夹角余弦 cos_theta np.dot(v1, v2) / (np.linalg.norm(v1)*np.linalg.norm(v2) 1e-12) ortho_err[j-1, i-1] abs(90 - np.degrees(np.arccos(np.clip(cos_theta, -1, 1)))) print(最大正交偏差(度):, ortho_err.max()) print(超5度网格占比(%):, (ortho_err 5).mean()*100)逻辑说明正交性检查的核心是切向量的点积点积越接近 0说明两条网格线夹角越接近 90 度。这个脚本在投影坐标系下用平面距离近似网格距在公里级时精度足够如果网格覆盖纬度范围很大应在投影平面坐标上计算不要在经纬度原始坐标里硬算。参数说明阈值取 5 度是我个人的工程习惯有些精细模拟项目要求 2 度以内。如果超标网格集中在局部常见做法是调整岸线边界条件重新生成网格如果大面积超标说明区域形状不适合用单块正交网格要考虑把模型域切分成子区域或者改用非正交骨架的模型。水深平滑通常用限幅处理常见做法是把相邻网格深度的比值限制在一个范围内例如新水深 相邻水深按差商限幅后再迭代几次。POM 自带的平滑工具或第三方前处理脚本都能做关键是不能过度平滑以至于把航道和沙脊抹掉。平滑后立刻重新做一次单一网格的静水压力检验用一个均匀温盐场、无强迫的环境跑 24 小时看流速是否在数值噪声水平这一步能提前暴露地形梯度引发的虚假流动。3.3 σ 分层设置层数、拉伸系数与浅水区厚度POM 的分层参数一般写在输入文件里常见需要指定总层数、垂向拉伸系数。总层数在工程实践里常见范围是 10 到 30 层近岸工程项目用 16 到 21 层较为常见。层数太少近底摩擦和温盐跃层分辨率不够层数太多外模虽然不受影响内模计算量近线性增加而且 σ 层在深水区会被拉得很厚对深层流动的解析未必有改善。拉伸系数决定层在表层和底部的加密程度。POM 的 σ 层默认是表底层加密、中间层较疏的分布这符合海洋中大部分流速和温盐变化集中在上层和底边界层的物理事实。参数表如下参数常见取值说明总层数16~21浅水工程模拟常用深水区可加到 26表层加密系数0.2~0.4越大表层越密风强迫响应更准但步长受限底层加密系数0.2~0.3决定底摩擦和底边界层解析能力最小水深 cutoff4~10 m低于此值按此值处理防浅滩层厚趋零底摩擦系数0.002~0.003量级参考需按底质类型率定逻辑说明σ 坐标的层厚是总水深除以层数浅水区一旦水深只有 5 米而层数 21表层厚度可能只有几十厘米垂向流速梯度大数值稳定性会变差。所以 POM 工程里普遍设置最小水深 cutoff把低于阈值的节点统一抬升到一个可计算深度。cutoff 设太高会改变岸线位置和潮间带面积设太低则层厚过小需要与底摩擦系数配合调整。参数说明层数改变后初始温盐场的垂直插值必须重新做不能沿用上一版的分层插值结果。拉伸系数看起来只是“密度分配”实际上它影响 Mellor-Yamada 湍流闭合在表底层的分辨率改完拉伸系数后要重新跑一遍静水检验确认没有在层与层之间产生虚假的垂向混合。4. 初始场、开边界与强迫场让冷启动模型不翻车的三件事4.1 初始温盐场冷启动还是气候态插值POM 是自由面原始方程模型对初始场不敏感是它相对一些全球模型的优点但“不敏感”不等于“可以随便给”。均匀温盐场的冷启动在深水区会产生明显的地转调整过程初始几天流场像噪声纯粹是模式在自我调整。如果这个过程发生在你关心的模拟窗口内会污染统计结果。常见做法是首月采用“气候态插值 一个星期的泄放期”把泄放期的结果丢弃。气候态温盐数据从公开的再分析或气候态数据集中获取插值到模型网格的 σ 层上。插值示意如下# 读取气候态温盐经度、纬度、深度、温度、盐度 ds nc.Dataset(clim_temp_sal.nc) depth_clim ds.variables[depth][:] temp_clim ds.variables[temp][:] # shape: (depth, lat, lon) salt_clim ds.variables[salt][:] # 定义模型 σ 层对应的实际深度以平均水深近似冷启动阶段足够 # sigma 层坐标通常为 -1 到 0深度 sigma * H sigma_vals np.linspace(-1, 0, 21) # 21 层示意 depth_model sigma_vals * dep_model.mean() # 对每个水平网格点做垂向一维插值 # 这里用 np.interp 沿深度方向逐点操作示意核心逻辑 def interp_vert(temp_profile, depth_clim, z_target): return np.interp(z_target, depth_clim, temp_profile) # 实际项目建议用 netCDF 分块循环避免一次性载入大数组逻辑说明σ 层深度随水深变化所以严格做法是每个水平网格点上按实际水深 Hη 计算各层深度再对温盐剖面插值。冷启动阶段用平均水深近似影响不大但正式模拟应逐点计算。注意深水区气候态数据可能到不了海底超出数据范围的深度要做“底部分层保持”而不是外推否则温盐剖面上会出现向下的虚假梯度。参数说明插值后还要做一步重要性检查——检查密度场是否稳定。如果气候态数据里的温度盐度经过插值后出现密度倒置模型会在前几步产生猛烈对流失稳。常见做法是插值后做静态稳定性修正把不稳定剖面的密度梯度调平然后重新计算温盐。4.2 开边界与潮汐驱动水位和流速怎么逼进模型POM 的边界处理常见组合是“辐射边界 边界强迫”。辐射条件的本意是让模型内部产生的波动能自由传出边界不被虚假反射。潮汐驱动的常规做法是把开边界节点上的水位按分潮叠加作为强迫输入常用分潮包括 M2、S2、N2、K1、O1振幅和迟角取自全球潮汐模型或区域潮汐数据集。实际操作中要设置两类边界信息外模边界每个开边界节点的水位时间序列或直接给调和常数由模型内部合成水位。内模边界温度、盐度、流速的边界值通常用气候态值加边界松弛方案。边界强迫的常见实现方式是把调和常数写成边界文件运行中每个外模步长合成一次水位。以 M2 分潮为例某个边界节点的水位是h(t) A_m2 * cos(omega_m2 * t - phi_m2) A_s2 * cos(omega_s2 * t - phi_s2) ...这个合成逻辑可以用 Python 脚本预生成水位时间序列也可以写成 Fortran 子程序由模型在线计算。工程实践中我更倾向在线合成因为边界文件如果按外模步长输出文件大且容易被磁盘 IO 拖慢。开边界位置的选择比边界条件本身更重要。常见经验是边界离关注区域至少留 10 到 20 个网格距的过渡带否则边界强迫的误差会直接进入研究区域。如果模型域里有多个开边界要给所有开边界一致的分潮常数否则会在域内人为制造潮汐不连续的虚假信号。4.3 强迫加载 ramp 与运行从冷启动到稳定推进风场、潮汐、河流入流全部在第一个时间步突然加载等于给模型一个阶跃激发会在开边界和岸边产生虚假的重力波。常见做法是给所有强迫乘一个 ramp 因子ramp min(1.0, t / T_ramp)T_ramp 常见 12 到 72 小时风场取短一些潮汐取长一些。! 示意在主要时间循环内施加 ramp ramp_factor min(1.0, time / 43200.0) ! 12 小时 ramp wind_stress wind_stress_input * ramp_factor eta_bc eta_tide * ramp_factor逻辑说明ramp 不是物理过程纯粹是数值启动手段。它让强迫从零逐渐加满避免冷启动阶段的高频振荡。ramp 结束后的前几个潮周期仍属于泄放期这一段的结果不要直接用于统计验证。参数说明ramp 时长没有通用最优值与模型域尺度、水深有关。小海湾 6 小时可能就够陆架尺度模型建议 48 小时以上。判断指标是自由面高度的振荡幅度如果 ramp 结束后水位序列里还能看到明显的高频波动说明 ramp 太短或外模步长不够安全。运行阶段看输出日志./pom pom.in run.log 21 # 实时观察推进是否正常 grep time run.log | tail -5 # 检查是否有 NaN grep -i nan run.log | head -5逻辑说明POM 运行日志会周期性打印当前时间步和部分诊断量。看到 NaN 不要急着改物理参数先回到时间步长和网格检查。日志里出现剧烈水位跳变也需要立刻终止跑批典型原因是外模步长过大、开边界强迫阶跃或者水深场有异常尖点。5. POM 建模避坑记录发散、振荡、正交性差的 4 个典型现场5.1 症状第一条时间步就 NaN日志里外模水位爆掉现象运行日志在第二十几个外模步长出现 NaN之前几个步长的水位数值交替出现巨大正负跳变。原因外模步长超出 CFL 条件。POM 外模求解的是二维浅水波波速约 √(gH)在深水区或深航道处波速快同样的空间步长下允许的时间步长更小。很多人只调了内模步长忘了外模步长也要同步缩放。解决先按最小网格距和最大水深估算步长上限再留安全系数。估算公式为 Δt Δx / √(gH_max)实际使用再除以 2 到 3。比如 Δx 500 m、H_max 60 m粗算上限约 20.6 秒工程上可取外模步长 5 到 10 秒内模步长对应取 50 到 100 秒。改完重新看日志。# 按最大水深估算外模步长上限并给出建议值 g 9.81 dx 500.0 # 最小网格距单位 m h_max 60.0 # 域内最大水深单位 m dt_upper dx / np.sqrt(g * h_max) dt_est dt_upper / 2.5 # 安全系数 2.5工程常用 print(fCFL 上界: {dt_upper:.1f} s, 建议外模步长: {dt_est:.1f} s)逻辑说明外模步长受限的本质是波动不能在一个步长内跨过多个网格。安全系数取 2.5 是我个人习惯有的项目取 3网格质量差、地形陡的区域建议取更大。这个估算只是初值跑起来后要看日志里是否有水位振荡。5.2 症状网格看着没问题湾内流场持续振荡现象静水检验就不过关无强迫环境下湾内出现几厘米每秒的持续流场形态固定不衰减。原因水深梯度过大引发 σ 坐标压力梯度误差以及网格正交性偏差超标。σ 坐标在陡峭地形上把水平压力梯度拆成两个大数相减计算机精度有限就产生了残留虚假力。误差大小正比于地形梯度与层数乘积。解决对水深场做限幅平滑逐网格限制相邻水深变化比在阈值 0.2 到 0.3 之间迭代平滑几轮同时检查正交性角度把超标网格的岸线控制点重新布置。两个动作做完重跑静水检验流速应降到毫米每秒量级。准则静水检验不通过就不要继续往下跑这是 POM 建模开始阶段的底线检查。任何后续结果的解释都必须建立在这一步干净的基础上。5.3 症状冷启动前 10 天结果像噪声统计指标全乱现象水位和流速时间序列前 10 天高频振荡明显第 10 天之后逐渐正常但前段数据又不能直接丢掉因为模拟窗口总共只有 30 天。原因强迫场阶跃加载风场和潮汐在第一个步长就施加全值冷启动状态下的密度场与强迫不匹配模式以重力波形式释放初始调整能。这不是模型 bug是启动方式问题。解决给所有强迫加 ramp潮汐 ramp 取 24 到 48 小时风场 ramp 取 6 到 12 小时同时延长泄放期至至少 3 个潮周期。如果项目时间窗口紧张把初始温盐场从气候态插值而不是均匀场能显著缩短调整期。经验不要指望通过减小时间步长来治好这种振荡步长只能缓解数值不稳定解决不了初始场与强迫失配的物理问题。5.4 症状浅水区 σ 层厚度异常近底层流速出现尖峰现象岸线附近的浅水网格点流速在特定潮相位出现尖锐脉冲单点时间序列上像毛刺空间上集中在一两个网格。原因最小水深 cutoff 设置过高导致实际计算水深远大于真实水深潮间带面积被压缩浅水区层数虽多但层厚极小垂向流速剪切被放大或者 cutoff 设太低σ 层在浅水点厚度只有几厘米数值截断误差被底摩擦放大。解决把最小水深 cutoff 调整为与潮差量级匹配例如平均潮差 2 米的区域cutoff 取 5 米左右再检查该处的 σ 层厚度是否小于一个阈值调整总层数或拉伸系数。单点尖峰也可能是底摩擦系数与层厚不匹配底摩擦系数在浅水区可以按空间变化场给出不要全域用同一个常数。注意改完 cutoff 要重新做岸线位置检查cutoff 抬高会把原本可见的潮间滩地变成永久淹没区对湿地岸段的水动力结果影响很大。6. 验证结果的一手习惯调和分析、守恒检查与 run 卡片模型跑通不是终点验证才是让结果能拿去用的关卡。验证 POM 结果我习惯抓三件事。第一是潮位验证把模拟水位做调和分析把 M2、S2、K1、O1 的振幅和迟角拿出来与观测水位站的调和常数对比计算振幅误差与迟角绝对误差。水位对不上流速基本也不用看因为潮汐是近岸流场的主驱动。调和分析可以用现成工具做最小二乘拟合示意逻辑如下# 用最小二乘拟合主要分潮振幅与迟角 t np.arange(0, 30*24*3600, 3600) # 逐小时模拟时间 eta_model water_level_time_series # 模型水位 periods {M2: 12.42*3600, K1: 23.93*3600, O1: 25.82*3600} A np.stack([np.cos(2*np.pi*t/T) for T in periods.values()] [np.sin(2*np.pi*t/T) for T in periods.values()] [np.ones_like(t)], axis1) coef, *_ np.linalg.lstsq(A, eta_model, rcondNone) # 由 cos/sin 系数换算振幅与迟角逻辑说明拟合长度至少覆盖 15 天才能分离 M2 和 S2 两个频率接近的分潮。与观测对比时迟角误差比振幅误差更敏感很多模型振幅率定得很好、相位差出 30 度说明底摩擦或水深分布仍不准确。第二是守恒检查。POM 在水位流量交换为零的封闭环境里理论上应严格守恒水量。在开边界模型里检查目标是长期模拟时域平均水位不发生持续漂移。如果 30 天模拟的水位均值在持续下降或上升通常是开边界强迫不平衡或外模与内模耦合存在净通量误差。第三是我自己的一个工作习惯维护一张 run 卡片表把每次运行的关键参数记录成一行。卡片的必需列包括网格版本、水深平滑轮次、总层数、外模步长、内模步长、底摩擦系数、cutoff 水深、ramp 时长、强迫数据版本以及这次跑出来的 M2 振幅误差。没有这张卡片参数调了几天之后回头看之前的“好结果”复现不出来那才是建模过程里最消耗信心的事。调参在这个模型上很容易变成玄学而 run 卡片是唯一的后悔药。这个习惯是在某个模拟项目上血泪换来的当时连续两天调底摩擦系数每次都觉得“这次看起来不错”最后发现记混了网格版本前面几组结果全都作废。从那以后每个 run 都先写卡片再提交验证完结果再归档。现在把这个习惯推荐给你。希望你建立模型的第一天就把它用起来希望帮到你。本文还有配套的精品资源点击获取
)
