FVCOM风场驱动文件制作：mkwndfv插值原理、参数与常见坑

简介面向FVCOM海洋数值模拟研究人员和海洋科学专业学生这份资源提供了一个专门用于风场数据预处理的MATLAB脚本。FVCOM作为非结构网格海洋环流模型对输入风场有着严格的时间和空间格式要求该脚本正是为了解决原始气象风场与模型需求之间的衔接问题使用者可通过简单的参数设置将它嵌入自己的工作流。资源压缩包内仅包含一个m脚本文件整体大小约1KB属于高度聚焦的小工具尽管体量小脚本覆盖了从数据读取、坐标转换、空间插值到非结构化网格、时间重采样、边界赋值以及输出FVCOM标准ASCII风场文件的完整逻辑并预留了用户自定义数据源与网格接口。目前已有761人浏览学习适合需要快速生成风场驱动文件、或希望借鉴成熟代码精简自研工具的建模工程师。通过阅读该脚本用户不仅能减少手工处理风场数据的重复劳动还能理解FVCOM风场文件的格式组织与质量控制要点为后续开展潮汐、环流或输运模拟提供可靠的大气驱动力。1. 为什么 FVCOM 风场要单独做mkwndfv 在做的事跑 FVCOM 做风暴潮或溢油模拟时最容易忽略的一步是把气象风场变成模型能吃的驱动文件。很多刚接触的人直接把再分析风场的经纬度格点喂给 FVCOM结果运行到一半风应力全变成 NaN或者边界上出现一条一条的锯齿状流速带。这个标题里的 mkwndfv就是 FVCOM 工具链里专门负责把格点风场插值到非结构三角网格节点上的预处理程序它解决的问题是风场源文件NCEP/GFS 的规则经纬度网格或者 WRF 输出的公里级网格和 FVCOM 的节点坐标是两套完全不同的离散系统必须先把风速值映射到每个节点上再按模型的风应力公式换算。这篇文章面向的是已经在跑 FVCOM、或者正准备把风应力加进自己的算例的工程师和海洋环境建模人员。我会从 mkwndfv 的选型逻辑讲起给出最小可运行的命令和参数表再把插值半径、时间单位、陆地点掩膜这几个最常踩的坑列出来。它能直接解决的问题是让你在半小时内生成一份风场驱动文件并且知道结果对不对。2. 选型与原理mkwndfv 与其他风场工具的边界2.1 先分清风场工具mkwndfv 不是唯一方案FVCOM 的历史版本里出现过多个风场相关工具名字相近但分工不同。常见的有mk_wind、mkwnd、make_wind和标题里的mkwndfv。它们的共同点是读取外部风场数据输出 FVCOM 需要的节点风速文件但输入格式和插值策略有差异。比较新的 FVCOM 版本通常推荐用make_wind或直接通过 WW3/耦合器传入风场mkwndfv是早期但依然在老旧项目里大量留存的工具。如果你的课题是接手别人跑了多年的风暴潮模型大概率会在预处理目录下看到mkwndfv的源码和 dat 控制文件。我一般这样判断该用哪个如果你的风场源是 WRF output且网格是旋转极地投影优先走make_wind的 WRF 分支如果源文件是 GFS 或者 NCEP 再分析这种等经纬度网格mkwndfv足够且更简单。它读入的是一份文本控制文件不用编译额外的 Python 接口对长期批量跑算例来说少了依赖就是少了一份故障源。2.2 mkwndfv 的插值逻辑反距离权重和时间线性插值mkwndfv的核心插值方法是反距离权重IDW对每个 FVCOM 节点在源风场网格上搜索若干临近点按距离倒数加权平均得到该节点的风速。这听起来简单但它有三个细节决定了结果第一搜索半径的选取。mkwndfv的控制文件里通常有spd或rmax参数代表最大搜索半径单位是度或者对应源网格的格点数。半径太小近岸节点会搜不到足够的风源点产生空白半径太大会把不同天气系统下的风场平滑掉近岸急流轴会被抹平。常见设置为源网格水平分辨率的 2 到 3 倍比如 0.25 度分辨率的 GFS取半径 0.5 到 0.75 度。第二时间维的线性插值。再分析风场通常是 3 小时或 6 小时间隔而 FVCOM 的驱动需要逐时或更密的风应力。mkwndfv会在两个时间层之间做线性插值并在输出文件中给出插值后的时间序列。这个过程的控制参数是time_dflt它决定输出的时间步长是多少秒或小时。第三风速到风应力的换算。mkwndfv输出的是风速还是风应力取决于编译时宏和源文件里是否要求做拖曳系数换算。多数版本输出的是 10 米风速u10/v10让 FVCOM 运行时按 GOTM 或 LSU 拖曳系数方案自行换算。这一点务必在读完头文件后确认否则把风应力当风速输入相当于把力放大了一个数量级。2.3 输入数据组织的惯例一个 dat 文件驱动一批格点数据mkwndfv的封装方式是典型的 Fortran 工具风格你准备好一份*.dat控制文件里面写清格点风场的路径、变量名、输出路径和插值参数然后命令行执行mkwndfv wind.dat或者mkwndfv wind.dat your_grid.grd。网格文件通常用 FVCOM 的节点坐标文件grd 格式里面是每个节点的经度、纬度和水深值。mkwndfv的处理对象是时间上连续的净 CDF 风场而不是单个时刻。所以你的源数据最好已经按时间维拼接成一个文件或者控制文件里支持列表方式指定多个文件。常见做法是先用 Python 把 GFS 按时间段合并成一个wind_input.nc再写 dat 文件指向它。3. 用最小数据跑通 mkwndfv准备、命令与输出3.1 准备一份 mkwndfv 能读的风场 netCDFmkwndfv的各版本对 netCDF 变量名的要求不完全一致但最宽松的常见约定是变量名uwnd/vwnd维度是time、lat、lon单位分别为 m/s 和度。如果你的源文件是标准 GFS 命名u10/v10最好先用 Python 统一一下。下面这段脚本把 GFS 的u10、v10转成uwnd、vwnd同时把经度从 0~360 转到 -180~180避免 mkwndfv 在跨日期变更线区域找不到点。import xarray as xr ds xr.open_dataset(gfs.t00z.pgrb2.0p25.f000.nc) da ds.rename({u10: uwnd, v10: vwnd}) # 经度归一到 [-180, 180) lon da[lon].values lon ((lon 180) % 360) - 180 da da.assign_coords(lonlon) da da.sortby(da[lon]) ds_out da[[uwnd, vwnd]].to_dataset() ds_out.to_netcdf(wind_input.nc)这段逻辑里有两处要注意rename只改名字不改坐标sortby是关键如果源文件的经度逆序排列插值时极点会乱。wind_input.nc里time的单位建议用hours since 1990-01-01因为 Fortran 工具对基准时间解析能力弱用通用历元能少踩坑。3.2 dat 文件模板与最小运行命令创建一个文本文件fvcom_wind.dat内容如下NCEP Wind File: wind_input.nc Grid File: fvcom_grd.dat Output File: wind_fvcom.nc Start Time: 2017-08-20 00:00:00 End Time: 2017-08-20 06:00:00 Time Step: 3600.0 Search Radius: 0.75 Missing Value: -9999.0然后运行./mkwndfv fvcom_wind.dat运行后终端会打印每个时刻的插值完成情况如果没有任何NaN警告说明源网格和节点坐标对上了。参数说明Search Radius单位是度控制每个 FVCOM 节点搜周围多少个源格点。0.75 对应 0.25 度 GFS 的三倍距离也是我常用的起手值。Time Step输出风场的时间间隔秒。3600 表示逐小时输出如果有潮位耦合需求可以压缩到 900。Missing Value源数据里的缺测填值一般设为 -9999.0插值逻辑遇到该值会在权重计算时跳过它。3.3 验证输出文件的基本结构插值完成后用 Python 检查输出文件里每个节点的风场是否合理以及时间维度是否和期望一致。import xarray as xr ds xr.open_dataset(wind_fvcom.nc) print(ds[u10].shape) # (time, node) print(ds[node].shape) # (node, ) print(ds[u10].values.flatten()[::10]) # 抽样打印节点值 # 检查陆地点是否为 0 或 nan import numpy as np node_slice ds[u10].isel(time0).values print(np.nanmin(node_slice), np.nanmax(node_slice))输出里node维应等于 FVCOM 网格节点数。如果形状里出现lat、lon维说明 mkwndfv 版本较旧输出的是格点化风场此时需要改成用scrip插值而不是直接给 FVCOM 读取。这个差异很多人查很久看到形状就能立刻判断。4. 参数与常见坑搜索半径、时间单位与陆地点掩膜4.1 三个必调参数Search Radius、Time Step、Missing Value 的语义这三个参数里Search Radius对模拟结果影响最大它直接决定近岸风场的空间结构。参数作用常见值影响Search RadiusIDW 搜索半径度0.5 ~ 1.5太小产生空白太大抹平风速梯度Time Step输出时间间隔秒900 ~ 3600影响模型对风应力变化的响应Missing Value源数据缺测值-9999.0不设对会找不到有效数据整场 NaN我特别强调一下 Search Radius 的边界如果源风场分辨率是 0.25 度取 0.5 有点紧遇到台风这种强梯度系统会在眼壁附近产生不连续的风速台阶取 1.5 以上则让风场像被高斯滤波过防御潮峰值偏低。如果你要复现台风期间的增水先用 0.75 跑再对比站点观测风逐步下调到 0.6 左右。4.2 时间单位错位的典型症状与绕过方案mkwndfv的老版本在读取 netCDF 时间变量时只接受hours since 1900-01-01之类的绝对基准如果源数据用了seconds since 1970-01-01读出来的时间会偏移约 8 小时取决于本地时区。症状是输出风场的时间标签和实际时刻对不上模型里风暴潮峰值提前或滞后一个时段。解决方法是统一在 Python 里转成时间变量字符串并显式写出基准时间的日期和时间示例ds[time].encoding[units] hours since 1990-01-01 ds.to_netcdf(wind_input.nc)如果输出的时间仍然错位检查你的 netCDF 库在写入时是否自动加了日历属性gregorian或proleptic_gregorian某些 Fortran 读入端会把默认日历当成360_day导致每月偏移 5~6 天。如果遇到这个情况可以在 Python 里删除calendar属性再写回。4.3 陆地点掩膜把插值后的风速归零FVCOM 自身的干湿网格可以在浅水区把陆地点自动处理掉但风场驱动文件里的陆地点如果仍保留风速值会在水位不高时产生虚假的风应力振荡。mkwndfv的多数版本会自动使用网格文件里的水深值判断水深为负或为 0 的节点风速置为 0但有的版本不做这个判断。我一般会在插值后用 NCO 工具或 Python 强制掩膜python mask_wind.py wind_fvcom.nc fvcom_grd.dat实现方法读 grd 文件的水深对水深 0 的节点把整个时间维的u10、v10置为零。注意不要使用NaNFVCOM 读入时对NaN的处理不稳定可能直接让整个扇区的计算挂掉。使用 0 是安全选择而且实际风应力公式在风速为零时输出为零不影响动量守恒。如果用 Python注意节点索引从 0 开始和 Fortran 的 1 开始差一。我曾经因为索引偏移把海里的风场抹掉一半排查了很久。4.4 排错速查三个最普遍的报警运行mkwndfv时出现FATAL: Node outside valid data range之类错误通常不是程序坏了而是源网格的经纬度范围没覆盖到 FVCOM 的边界。检查源风场的min(lat)是否小于你的研究区最南端max(lat)是否大于最北端。GFS 全球数据基本不会出这问题但 WRF 区域数据很容易出边界。第二常见的问题是输出风场出现条带状图案。原因是源风场网格在某些区域缺测IDW 搜索半径内找不到足够多有效点只能靠少数几个点强行插值。解决方法是把Missing Value检查打开并在 dat 文件里允许忽视缺测行。第三是时间维混乱常见原因是多个源文件拼接时时间不连续。正确做法是在 Python 里重新生成时间坐标然后对所有变量同步排序ds ds.sortby(ds[time])这个看起来基础但批量处理多年数据时经常因为glob顺序问题把 2018 年排在 2019 年后面导致插值出的时间序列突变跳变。5. 进阶多源风场拼接与逐时验证的 3 个脚本技巧5.1 用列表批量合成长期风场做 10 年以上模拟时你不会想要手动写 dat 文件指向几十个季度的数据。常见做法是写一个脚本生成fvcom_wind.dat的片段每次运行只改文件路径和时间范围。for year in 2018 2019 2020; do sed s/START/${year}-01-01/; s/END/${year}-12-31/ wind_template.dat wind_${year}.dat ./mkwndfv wind_${year}.dat donesed只是批量替换时间范围输出文件按年份命名最后用 CDO 或 Python 把多个wind_fvcom_.nc合并。注意合并时要做时间排序和去重否则 FVCOM 的读入端会报时间不单调。5.2 验证插值半径对岸线风应力的影响要量化不同半径对模拟结果的影响可以在相同的日期范围内跑三组半径然后对比研究区内某个站点附近的节点风速时间序列。使用xarray先取离站点最近的节点import xarray as xr import numpy as np site_lon, site_lat 122.0, 30.5 for r in [r05, r075, r10]: ds xr.open_dataset(fwind_fvcom_{r}.nc) dist (ds[lon] - site_lon)**2 (ds[lat] - site_lat)**2 idx int(dist.argmin()) u ds[u10][:, idx].values np.savetxt(fu10_{r}.txt, u)比较三份文本文件的方差和峰值时刻如果峰值差异超过 15%就说明源网格太粗糙此时应该改用更高分辨率风场或调整插值方案而不是继续调半径。5.3 用潮位站点风记录做快速交叉验证在投入 FVCOM 运行前把插值出的节点风场与海岸气象站观测做散点对比是性价比最高的验证。步骤如下从输出文件里提取离站点最近的节点风速。把观测站的 10 分钟平均风换算成 10 米高度、逐时平均。计算平均绝对误差和相关系数。对风速而言相关系数达到 0.85 以上说明插值合理如果低于 0.7优先检查站点附近的地形是否未在网格中体现其次检查源风场的空间分辨率。这个验证比跑到最后看增水误差容易得多也能在早期定位风场文件本身的问题。本文还有配套的精品资源点击获取