简介面向FVCOM海洋模型用户的MATLAB风场预处理脚本解决为FVCOM非结构化网格准备风场输入数据的问题。脚本涵盖读取原始风场、空间插值、时间重采样、坐标转换到输出FVCOM可识别ASCII文件的完整处理流程支持按需调整物理参数与文件命名规则能够适配GFS、ECMWF等再分析风场产品也可扩展处理站点观测数据。压缩包仅含1个m文件大小约1KB轻量易用适合海洋、气象领域研究人员及FVCOM初学者参考。已有761人学习下载。通过该脚本可深入理解FVCOM风场制备的关键环节包括数据获取与质量控制、从规则网格到非结构化网格的空间插值、模拟时间步长匹配、开边界与陆边界条件赋值等同时可作为模板二次开发快速迁移到其他海域或不同风场数据源显著减少手工编写数据接口的时间成本为后续FVCOM模拟提供规范可靠的气象驱动输入。1. 风场数据准备为什么总卡在 mkwndfv 这一步跑 FVCOM 的人大概率都经历过同一个尴尬模式装好了、网格画好了、初始场也转出来了结果卡在大气强迫场这一步。风场数据格式不对、节点顺序对不上、时间步长不匹配任何一个问题都能让模式直接跳出 NaN。mkwndfv.m这个脚本解决的就是这最后一公里——把 GFS、ERA5 这类再分析风场转成 FVCOM 能直接读取的三维风速强迫文件。适合的读者不是写模式的人而是被数据预处理耗掉两三天、只想知道“到底哪一步错了”的人。这篇文章会带你拆开mkwndfv.m的完整数据流从 NetCDF 里抽出 U/V 分量通过插值落到 FVCOM 非结构化网格节点上再按时间步长输出成 ASCII 文件。中间穿插我实际跑过的参数配置和踩坑记录最后给出一套验证风场质量的脚本思路。2. 读懂 mkwndfv.m 的输入输出从再分析风场到 FVCOM ASCII2.1 风场源数据选型GFS、ERA5 与 FVCOM 网格的差异FVCOM 本身不对风场做任何插值它只负责在每步计算时去读你给的文件。这意味着你给它什么分辨率、什么投影、什么变量名它就认为那已经是模型网格节点上的值。常见的风场源有三种它们的特性和适配场景有明显区别。GFS 是 NCEP 的全球预报产品0.25 度分辨率每 3 小时或 6 小时一个时次适合区域较大、对时效有要求的业务化模拟。ERA5 是 ECMWF 的再分析资料0.25 度逐小时输出质量更稳但文件体量也大下载和读取都更耗时。如果只做局部海域的高分辨率模拟还可以用 CFSR 或本地的 WRF 输出此时坐标系通常已经是 Lambert 或 Mercator需要额外做投影转换。选择时先看你的 FVCOM 网格范围。如果模拟域只有数百公里用 GFS 的原始网格直接插值会产生明显的边界锯齿建议先做一个中间网格的平滑。另外注意风场源文件的变量名GFS 里是 UGRD/VGRDERA5 是 u10/v10读数据前用ncinfo确认否则索引取错后面所有的插值都是白算。2.2 mkwndfv.m 的输入参数与调用方式mkwndfv.m在压缩包里是独立文件不依赖额外的工具箱函数基础 MATLAB 环境就能跑。它典型的调用方式是mkwndfv(wind_era5_202306.nc, fvcom_mesh.dat, outdir, 6)第一个参数是原始风场 NetCDF 文件路径第二个是 FVCOM 网格文件路径outdir指定输出目录最后的6表示输出时间间隔小时。如果原始数据是逐小时的而模式步长要求每 6 小时一个风场文件脚本内部会做均值或最近时次抽取。参数含义并不复杂但有几个细节值得注意。第一fvcom_mesh.dat是 FVCOM 的标准网格文件里面包含节点数、单元数、节点坐标和深度。脚本通过load或文本解析读取节点坐标而不是从 NetCDF 里读所以你必须保证这个文件和模式运行时用的是同一个版本。第二原始风场的经度范围如果是 0–360而 FVCOM 网格是 -180–180脚本里必须有经度平移的判断否则插值结果是全零。我一般会在调用前先检查网格文件的坐标单位FVCOM 常用的是经纬度十进制度但也有用 UTM 的情况。如果网格是 UTM而风场是经纬度需要先把风场坐标投影到 UTM再插值到节点。很多版本没有内置这个功能需要你提前用m_map或者deg2utm转换好坐标再传给脚本。2.3 输出文件规范节点编号、U/V 分量与时间文件名FVCOM 的大气强迫风场通常是 ASCII 格式每个时间步对应一个独立文件文件命名直接决定模式能否按顺序读取。例如wind_20230601_00.dat、wind_20230601_06.dat。文件内容形如5181 1 -3.42 5.18 2 -3.01 4.77 ... 5181 -2.88 4.31首行是节点数后面每一行依次为节点编号、U 分量东西向m/s、V 分量南北向m/s。注意这里的 U 和 V 是地球坐标下的不是旋转到 FVCOM 网格切向的法向分量。部分版本要求输出风速和风向两个字段此时第三列是风速第四列是风向度。你的mkwndfv.m用的是哪种格式直接决定 NML 文件里的读取方式输出前最好打开一个文件确认。输出路径和文件名里的时间是 UTC不是本地时间这一点容易被忽略。很多跑业务预报的人习惯用北京时导致风场提前或滞后了 8 小时。处理方式是文件命名统一用 UTC而真正算模式的起始时间也在 NML 里写 UTC这样不管你在哪个时区都按模式内部时间走。3. 核心代码拆解插值、旋转与时间重采样3.1 读取 NetCDF 风场并提取 U/V 分量mkwndfv.m的核心逻辑并不复杂难的是每一步都要考虑数据源的差异。读取 NetCDF 我常用ncread先通过ncinfo拿到变量名和维度信息再动态提取。ncfile wind_era5_202306.nc; u10 ncread(ncfile, u10); v10 ncread(ncfile, v10); lon ncread(ncfile, longitude); lat ncread(ncfile, latitude); time ncread(ncfile, time);这段代码里u10和v10的维度一般是[lon, lat, time]注意 NetCDF 的维度顺序可能不同需要用ncinfo确认。ERA5 的长期变量是longitude和latitudeGFS 里则可能是lon和lat直接写死会报错所以我喜欢加一层判断if isfield(info, lon) lon ncread(ncfile, lon); else lon ncread(ncfile, longitude); end这样至少能兼容大多数常见数据源。时间变量time的单位是hours since 1900-01-01之类的需要换算成datenum才能和 FVCOM 的模拟时间对齐。3.2 从 FVCOM 网格文件读取节点坐标FVCOM 的网格文件也叫casename_grd.dat格式相对固定。首行是节点数然后每个节点一行包含节点编号、x 坐标、y 坐标、水深、类型标记。读取时用文本扫描比load更稳因为还混着单元数信息。fid fopen(fvcom_mesh.dat, r); nnode str2double(fgetl(fid)); grid_data textscan(fid, %d %f %f %f %d, nnode); node_id grid_data{1}; x_node grid_data{2}; y_node grid_data{3}; fclose(fid);注意textscan的格式字符串要和实际文件列数严格对应否则后面的单元部分会被误读。如果网格文件里每行还有额外的标记位比如边界标记应改为%d %f %f %f %d %d来吸收多余列。读取后最好检查x_node的数值范围判断是经纬度还是米制坐标。3.3 双线性插值与掩膜处理风场源数据是规则网格而 FVCOM 节点是非结构分布因此核心步骤是把规则网格上的值插到每个节点上。MATLAB 自带的griddata可以用但速度偏慢而且边界处容易产生外插值。我一般优先用interp2因为它更快、更可控。[X, Y] meshgrid(lon, lat); U_interp interp2(X, Y, u10(:,:,1), x_node, y_node, linear);这里有个隐藏问题源风场的经纬度矩阵如果是[lat, lon]排列u10(:,:,1)的行对应纬度列对应经度interp2要求 X 和 Y 是单调递增的网格坐标所以必须把u10转置后再赋给X和Y。否则结果会出现明显的条带状噪声。另外interp2对超出源网格范围的点返回 NaN需要用掩膜把陆地节点或者外部节点的风速强制置零。陆地掩膜有两种做法。第一种是利用 FVCOM 网格文件里的水深值水深为负且绝对值小于 0 的节点视为陆地风速直接置零。第二种是使用源数据自带的lsmland-sea mask变量将海陆掩膜插值到节点上再对陆地节点赋值。推荐用第二种因为网格文件里的水深和风场数据的海岸线不是同一套坐标容易出现近岸节点漂移。3.4 输出 ASCLL 文件并检查字段顺序插值完成后输出阶段要把 U、V 按节点编号顺序写进文件。如果节点编号不连续或者源文件里是从 0 开始都要重新编号。输出时我一般写成定长格式省空间也方便读取。filename sprintf(wind_%s_%02d.dat, datestr(time_out, yyyymmdd), hour_out); fid fopen(fullfile(outdir, filename), w); fprintf(fid, %d\n, nnode); for i 1:nnode fprintf(fid, %d %8.4f %8.4f\n, node_id(i), U_node(i), V_node(i)); end fclose(fid);这里time_out是重采样后的时间hour_out是当天的时次比如 0、6、12、18。输出的 ASCII 文件不要包含任何注释行FVCOM 只认数字和空格多一个字符都会导致读取错位。写完文件后可以立即用fgetl读回来验证节点数避免到模式里才发现行数不对。4. 时间对齐与边界条件让风场和模式步长不打架4.1 时间维度重采样的两种思路FVCOM 的NML文件里WindForce部分会指定风场文件的频率通常和外部强迫时间间隔一致。如果你的源风场是 1 小时间隔而模式设定每 3 小时读一次风场那么有两种选择抽取整点小时还是做时间平均。时间抽取的优点是简单直接按索引选取即可但会丢失中间变化。时间平均更平滑不过会引入相位滞后。对于潮汐主导的河口区域风场的高频变化影响较小时间平均问题不大如果是台风过境场景峰值风速被平均后会削弱导致增水被低估。我一般建议保持源数据原始频率通过修改 NML 的WindFileInterval来匹配而不是牺牲数据质量。下面是一个简单的时间重采样代码按小时索引抽取最近的时次target_times datenum(2023,6,1,0:6:24,0,0); [~, idx] min(abs(time_num - target_times(1)));这里time_num是 NetCDF 时间变量转成的datenum序列target_times是你要输出的时间点。min(abs(...))找最近时次适合源数据时间不规则的情况。如果源数据本身是规则间隔直接用idx (target_times - time_num(1))/(time_num(2)-time_num(1)) 1更快。4.2 陆地边界与开放边界的风速处理FVCOM 模拟域通常包含陆地边界和开放海洋边界。风场在陆地上没有物理意义但你插值到 FVCOM 节点时如果节点位于陆地会出现风速不为零的荒唐结果。正确的处理方式是在插值完成后对所有陆上节点的风速置零。开放边界则不需要特殊处理因为风场本来就是全域的边界节点只要有值就行。但要注意一点如果 FVCOM 模型本身有干湿网格wet/dry那风场文件里的陆地节点需要设为 0否则干网格在蒸发和风应力计算时会出现异常。这里给出一个掩膜应用片段mask interp2(X, Y, lsm, x_node, y_node, nearest); U_node(mask 0.5) 0; V_node(mask 0.5) 0;lsm建议用 0/1 二值量1 表示海洋。interp2用nearest避免在海陆交界处产生平滑过渡因为平滑值会导致近岸节点风速偏低进而影响风应力计算。4.3 时间文件命名和 NML 配置不一致的排查很多次模式报错不是风场内容不对而是文件名匹配不上。FVCOM 对时间文件的命名有严格约定例如wind_20230601_00.dat和 NML 里WindFileName wind_后缀自动拼接。如果你的脚本输出的文件名是wind_20230601_0.dat少了一个零FVCOM 就读不到。排查顺序是先看 NML 里WindFilePrefix和WindFileInterval再看实际文件名最后看文件里的节点数和网格文件是否一致。一个快速检查命令head -3 wind_20230601_00.dat wc -l wind_20230601_00.datwc -l的行数应该等于节点数加 1首行节点数如果多出来说明文件里有空行或注释行需要清理。如果少一行说明最后一个节点没写进去检查循环边界是否nnode而不是nnode-1。5. 一小时搞定风场验证连续性与物理合理性检查风场文件生成后不要急着跑模式先用脚本检查三个层面数值分布是否合理、空间场是否平滑、与观测对比误差能不能接受。我把这套验证写成了一个独立的 MATLAB 脚本每次生成风场后直接运行输出全是可视化图肉眼扫一遍就能发现数据源或插值的错误。5.1 空间分布检查一眼看出插值错位nodelon x_node; nodelat y_node; scatter(nodelon, nodelat, 6, U_node, filled); colorbar; title(U component at nodes); axis equal; xlabel(Longitude); ylabel(Latitude);运行后如果出现颜色分区明显错位根源大概率是经纬度排列顺序或者网格坐标单位没有统一。还有一种情况是边缘出现蓝色大块说明interp2的外插NaN被置零需要检查掩膜范围是否覆盖了模拟域。5.2 与站点观测对比的物理合理性假设你有浮标站点的实测风速数据可以提取离站点最近节点的模式风场值做对比。关键指标是均方根误差RMSE和偏差连续几天的对比能暴露系统性偏低或偏高的问题。obs_time datenum(2023,6,1:5,0,0,0); node_idx knnsearch([x_node y_node], [site_lon site_lat]); model_wind squeeze(U_node(node_idx,:)); rmse_val sqrt(mean((model_wind - obs_wind).^2));knnsearch找最近节点是常用做法注意节点坐标要用同一投影下的量。如果 RMSE 大于 3 m/s先检查是否拿 GFS 的 10m 风直接和站点 2m 风对比高度不同本身就有差异一般需要按对数风廓线换成 10m 标准高度。5.3 批量验证脚本一次跑完全部时间文件与其每个文件单独打开看不如写一个循环把所有输出文件读完绘制时间序列且统计缺失值。下面这段脚本按文件名顺序读取并检查有没有全部是零或连续 NaN 的时刻files dir(fullfile(outdir, wind_*.dat)); for k 1:length(files) fid fopen(fullfile(outdir, files(k).name), r); n str2double(fgetl(fid)); data textscan(fid, %d %f %f, n); fclose(fid); u_all data{2}; v_all data{3}; if all(u_all 0) all(v_all 0) fprintf(Time %s all zero!\n, files(k).name); end if any(isnan(u_all)) fprintf(Time %s contains NaN\n, files(k).name); end end这个循环跑完最快十几秒能立刻定位到异常时次。如果某个时刻全是零多半是源数据在该时次有缺测ncread返回了 fillvalue需要回到原始文件检查。如果只有近岸节点为零而外海正常那说明掩膜范围偏大把海洋节点也抹掉了。建议把这段脚本保存为check_wind.m以后每个新风场目录都丢进去跑一遍再决定是否启动 FVCOM 运算省去模式报错再回头查的时间。本文还有配套的精品资源点击获取