简介POM普林斯顿海洋模型是一套面向海洋数值模拟初学者与科研人员的开源建模工具可在f77、f90、g77、pgf77、CVF等多种编译器下直接运行适合用于河道水流、河口及近岸海域的水动力与温盐场模拟研究。资源包共21个文件约658KB包含f与c源码、n与dat输入数据、out运行输出、runpom2k执行脚本、change版本记录、pdf中文说明书及zip/tar工具包等覆盖从编译、算例运行到结果查看的完整流程。POM2K提供两个可直接上手的实例河道模型参数已内置于代码中河口模型则附带原始网格文件其输入文件由GRID-DATA/GRID.f生成便于用户理解网格与初始场的构造方式。借助GRID-DATA中的插值代码可将地形、温盐、表面风场及热通量等资料插值到网格点生成初始条件与强迫条件为自建模型提供有效工具。目前已有997人学习下载适合希望快速搭建POM建模环境、掌握算例配置与数据前处理思路的读者参考。1. POM 普林斯顿海洋模型从零建立一套可跑通的海域模拟如果你手头只有一份岸线数据、几组风场和开边界水位却要在一周内搭出一套能出流场和温盐结构的近岸模型POMPrinceton Ocean Model大概率是你绕不开的选择。它属于三维原始方程海洋模型采用 sigma 坐标贴合地形垂向用二阶湍流闭合时间积分走的是外模式与内模式分离的显式差分。这套结构决定了它对近岸、河口、陆架海域的适应性很强也决定了它对时间步长和网格质量相当敏感。我见过不少人拿到源码包后卡在编译和参数配置上其实真正难的不是代码本身而是搞清楚每个输入文件对应哪块物理过程。这篇笔记按“资源是什么、怎么配、坑在哪”的顺序拆一遍适合做近岸水动力、温盐输运、工程海域论证的从业者也适合刚接触海洋数值模拟的研究生照着复现。2. POM 的坐标系与差分框架为什么它适合近岸而不是深海2.1 sigma 坐标与内外模式分离的取舍POM 最核心的设计是垂向 sigma 坐标也就是把水深归一化到 0 到 1 之间自由表面和底部地形都能被贴合。这样做的好处是近岸浅水区的垂向分辨率不会被水深变化吃掉河口和陆架坡折带都能保持同样的垂向层数。代价是斜压梯度项在陡坡处容易产生误差所以实际配置时地形坡度一般控制在 0.05 以内超过这个值要么平滑地形要么加密水平网格。时间积分上POM 把正压模态和斜压模态拆开。外模式负责自由表面和垂向平均流时间步长受 CFL 条件限制通常取几秒到十几秒内模式负责温盐和三维流场步长可以放大到外模式的几十倍。这种分离让计算量降下来但也带来一个隐患内外模式之间的耦合频率如果设得不合理会出现质量不守恒或者表面重力波失真。常见做法是让内外模式步长比控制在 20 到 50 之间具体看网格分辨率和最大水深。2.2 控制方程与湍流闭合的落地含义POM 求解的是 Boussinesq 近似下的三维原始方程包括动量方程、连续方程、温度盐度输运方程和状态方程。动量方程里保留了科氏力、压强梯度力、湍流黏性项和底摩擦。湍流闭合用的是 Mellor-Yamada 二阶闭合方案也就是常说的 2.5 阶模型它通过求解湍动能和湍流宏观尺度两个预报量来给出垂向混合系数。对使用者来说这套闭合方案意味着你不需要自己调混合层深度但需要关注两个参数一个是湍动能产生项里的常数另一个是壁面函数里的粗糙度。前者在源码里通常以常量形式出现不建议随意改后者和底摩擦系数直接相关近岸模拟里一般取 0.002 到 0.005 之间。如果模拟出来的底层流速偏大先查底摩擦再查垂向分层是否被过度混合。2.3 网格与地形文件的准备步骤POM 的输入网格是正交曲线网格水平方向用 Arakawa C 网格离散也就是速度分量分别位于网格边中点标量位于网格中心。准备网格时第一步是确定研究区域和开边界位置第二步是生成正交曲线网格第三步是把水深插值到网格点上并做平滑。常见做法是用正交曲线网格生成工具先出网格文件再导出 POM 需要的 fortran 格式。水深数据一般来自实测或公开海图插值后要做三次平滑每次平滑系数不超过 0.5否则会抹掉真实地形特征。下面是一个把水深数据插值到 POM 网格并平滑的 Python 片段用的是常见科学计算库。import numpy as np from scipy.interpolate import griddata # 读取原始水深散点lon, lat, depth raw np.loadtxt(depth_raw.txt) lon, lat, depth raw[:, 0], raw[:, 1], raw[:, 2] # 读取 POM 网格点坐标 grid np.loadtxt(pom_grid.txt) gx, gy grid[:, 0], grid[:, 1] # 插值到网格点近岸用线性外海可换 cubic depth_grid griddata((lon, lat), depth, (gx, gy), methodlinear) # 填补插值空洞避免陆地边界出现 NaN mask np.isnan(depth_grid) if mask.any(): depth_grid[mask] griddata( (lon, lat), depth, (gx[mask], gy[mask]), methodnearest ) # 三次平滑系数 0.3保留地形趋势 smoothed depth_grid.copy() for _ in range(3): smoothed[1:-1, 1:-1] ( smoothed[1:-1, 1:-1] * 0.4 (smoothed[:-2, 1:-1] smoothed[2:, 1:-1] smoothed[1:-1, :-2] smoothed[1:-1, 2:]) * 0.15 ) np.savetxt(depth_pom.txt, smoothed, fmt%.3f)这段代码的逻辑是先做散点插值再处理边界空洞最后做三次加权平滑。参数上平滑系数 0.3 是经验值地形复杂的河口可以降到 0.2开阔陆架可以升到 0.4。插值方法选 linear 是为了避免 cubic 在数据稀疏区产生过冲如果原始水深点足够密可以换成 cubic 提升光滑度。输出格式按 POM 读取要求写成纯文本每行一个网格点的水深值。3. 编译与输入文件配置把源码跑起来的关键几步3.1 编译环境与 Makefile 调整POM 源码通常是 Fortran 77 和 Fortran 90 混合编译前先确认编译器支持固定格式和自由格式混编。常见做法是用 gfortran在 Makefile 里指定-ffixed-line-length-132和-fno-second-underscore避免长行被截断和符号名不匹配。如果用的是 ifort对应选项是-extend-source和-assume nounderscore。编译顺序一般是先编公共块和参数模块再编主程序。源码包里通常有一个makefile或Makefile里面已经写好了依赖关系但路径需要按本机环境改。我一般会先跑一遍make clean再make看第一个报错出现在哪个文件。如果是找不到netcdf.inc说明输出用了 NetCDF需要先装好库并在 Makefile 里加-I和-L路径。# 清理旧目标文件 make clean # 指定编译器和 NetCDF 路径后编译 make FCgfortran \ FFLAGS-O2 -ffixed-line-length-132 -fno-second-underscore \ NETCDF_INC/usr/local/include \ NETCDF_LIB/usr/local/lib -lnetcdf这里-O2是常规优化调试阶段可以换成-g -O0方便定位数组越界。NETCDF_INC和NETCDF_LIB按实际安装路径改如果不用 NetCDF 输出可以把相关源文件从 Makefile 里注释掉。编译成功后会在当前目录生成可执行文件名字通常是pom.exe或pom2k。3.2 输入文件清单与参数对应关系POM 的输入文件一般包括网格文件、水深文件、初始温盐文件、开边界文件、风场和热通量文件。每个文件对应一个 Fortran 读取单元单元号在源码的params或infile里定义。下面这张表列出常见文件和关键参数方便对照检查。文件名内容关键参数常见取值grid.dat网格坐标网格点数、正交性由生成工具决定depth.dat水深最小水深不小于 2 minit_ts.dat初始温盐表层温度、盐度按季节实测obc.dat开边界水位、流速、温盐潮汐或实测wind.dat风场风速、风向10 m 高度heat.dat热通量短波、长波、感热按月份配置时最容易出错的是开边界。POM 支持辐射边界和指定边界两种辐射边界适合外海开边界指定边界适合有实测潮位的断面。如果开边界给的是潮汐调和常数需要在源码里加一个潮汐预报模块或者提前把水位时间序列算好写进obc.dat。风场和热通量的时间步长要和模型时间步长对齐否则会插值出错。3.3 时间步长与 CFL 条件的实际估算时间步长的选取直接决定模型能不能稳定跑完。外模式步长受 CFL 条件限制公式是dt dx / sqrt(g * H_max)其中dx是最小网格间距H_max是最大水深g是重力加速度。举个例子最小网格 500 米最大水深 100 米算出来dt要小于 15 秒左右。实际配置时取 10 秒留余量。内模式步长通常是外模式的 20 到 50 倍但不要超过 60 倍否则斜压调整跟不上。如果模拟区域有强潮汐内外步长比取小一点比如 20如果是陆架环流可以取到 40。跑之前先用理想地形做一次冷启动看自由表面是否发散再换真实地形。# 冷启动检查先用 10 秒外步长跑 1 天 ./pom.exe run_control.txt # 查看输出日志里的 CFL 数和表面水位极值 grep -E CFL|SSH pom.log | tail -20如果日志里 CFL 数持续大于 1说明步长偏大或者地形有尖峰需要减小步长或平滑地形。表面水位极值如果出现异常大值先查开边界是否给了不合理的水位再查底摩擦是否过小。4. 避坑与排查跑 POM 时最容易翻车的五个地方4.1 编译通过但运行报数组越界现象是编译没问题一运行就Segmentation fault或者Fortran runtime error: Index out of bounds。原因通常是网格文件里的点数比源码里声明的数组维度大或者水深文件的行数和网格点数不一致。解决方法是先检查params文件里的IM和JM是否和网格文件头一致再逐行核对水深文件的行数。我一般会写一个简单的 Python 脚本统计行数和网格文件对比。4.2 自由表面发散或水位异常现象是跑了几十步后表面水位指数增长日志里 CFL 数越来越大。原因可能是外模式步长超过 CFL 限制或者开边界水位给了不合理的突变。解决方法是先把外步长减半试跑如果稳定了说明是步长问题如果还发散检查开边界文件的时间序列是否有跳变必要时做线性插值平滑。4.3 温盐出现负值或超过物理范围现象是温盐场跑着跑着出现负盐度或者温度超过 40 度。原因通常是平流项用了中心差分在强梯度区产生数值振荡。POM 里可以选择迎风差分或者加水平扩散。解决方法是把平流方案改成迎风或者把水平扩散系数调大但不要超过 100 m²/s否则会过度抹平锋面。4.4 底摩擦系数取错导致底层流速失真现象是表层流速合理底层流速偏大或偏小。原因是底摩擦系数没有按底质类型调整。沙质海底取 0.002 到 0.003泥质海底取 0.004 到 0.005。如果模拟区域底质不均可以按区域给不同系数但要在源码里加判断。我一般会先用一个中间值跑一遍再和实测流速对比微调。4.5 输出文件太大导致磁盘写满现象是跑了一半程序卡住日志没有报错但输出目录不再更新。原因是 NetCDF 或二进制输出没有做时间降采样每个时间步都写。解决方法是在输出模块里加一个输出间隔比如每 24 个内模式步写一次或者只输出表层和底层。如果已经写满先清理旧文件再改输出频率重跑。5. 进阶技巧用理想实验验证配置再上真实海域5.1 用理想地形做冷启动验证真实海域配置复杂直接上容易把配置错误和物理问题混在一起。我习惯先用一个理想化的矩形海域做冷启动地形是平底或者简单斜坡开边界给一个正弦潮汐初始温盐是常数。跑一天后检查表面水位是否和开边界同相位流速是否对称。如果理想实验都跑不通真实海域更不用谈。import numpy as np import matplotlib.pyplot as plt # 读取理想实验的表面水位输出 ssh np.loadtxt(ideal_ssh.txt) time np.arange(len(ssh)) * 3600 # 假设每小时输出一次 # 画水位时间序列检查是否和开边界同相位 plt.plot(time / 3600, ssh, labelmodel SSH) plt.axhline(0, colorgray, linestyle--) plt.xlabel(Time (hour)) plt.ylabel(SSH (m)) plt.legend() plt.savefig(ssh_check.png, dpi150)这段代码用来快速检查水位相位。如果模型水位和开边界相位差超过半个周期说明开边界辐射条件有问题需要检查辐射速度参数。参数上辐射速度一般取当地重力波速浅水区约等于sqrt(g * H)。5.2 用实测数据做温盐剖面校验理想实验通过后换真实地形和实测开边界。跑完一个潮周期后把模拟的温盐剖面和实测 CTD 数据对比。常见做法是取同一位置同一时刻的剖面画在同一张图上。如果表层温度偏差超过 1 度先查热通量文件的时间插值如果底层盐度偏差超过 0.5查底摩擦和垂向混合系数。校验项可接受偏差优先排查表面水位0.1 m开边界、底摩擦表层温度1.0 度热通量、短波穿透底层盐度0.5垂向混合、底摩擦表层流速0.1 m/s风场、开边界5.3 长期运行的稳定性习惯跑长期模拟时我一般会先跑 10 天检查能量是否守恒再决定是否延长。能量检查的方法是看动能和势能的日变化是否合理如果总能量持续增长说明有数值不稳定。另外每次改配置后都从冷启动重跑不要热启动叠加否则错误会累积。从那以后我每次改完参数都强制走一遍理想实验确认没问题再上真实海域。希望帮到你。本文还有配套的精品资源点击获取