简介gcmfaces是面向Matlab与Octave用户的开源工具箱专门用于处理全球气候模型GCM中的海洋环流数据。它采用face分块结构管理大规模网格数据提供数据读取、二维三维可视化、物理量计算与并行计算等功能适合从事海洋科学与气候研究的科研人员、研究生及有数据分析需求的工程师。压缩包内含317个文件以296个m脚本为主另有RST说明文档、PDF指南、YAML配置文件及少量辅助脚本整体仅3MB目录结构清晰。目前已有118人关注学习。通过这份资源读者可获得完整的工具箱源码、示例诊断模块、常用数据读取函数以及相关文档能帮助快速搭建海洋环流数据处理环境掌握分块网格数据的分析与可视化方法并支持在个人电脑或集群环境中高效处理高分辨率模型输出显著提升科研效率。1. 为什么海洋环流数据需要 gcmfaces从 face 分块说起如果你处理过 MITgcm 的全球海洋输出大概率见过那些名为*.001.001、*.002.001的二进制文件。它们看起来像一堆互不相干的小矩阵实际上是一个六面球面网格的六个面。直接用 Matlab 的load或squeeze去读得到的是连续但边界错乱的数组画图时等值线会在 face 交界处断裂。gcmfaces 就是为解决这个问题设计的 Matlab/Octave 工具箱它用一套重载后的数据结构和运算符让 face 分块数据可以像普通 lon-lat 矩阵一样做切片、平均、通量计算和绘图。它不是通用计算库而是专门面向 GCM 海洋环流诊断的专用层。下面基于源码包里的make.bat、rdmds.m、diags_set_*.m、process2nctiles.m等文件从数据读取到批量诊断完整过一遍适合正在跑 MITgcm 或处理 ECCO 系列输出的科研人员和工程师。2. 拆开 gcmfaces 包目录、核心函数与数据流先说结论这个包不是一个独立的软件而是一个 Matlab/Octave 工具箱核心是gcmfaces对象类和对 face 结构的重载运算。下载解压后你会发现根目录下直接放着make.bat、bibli.bib、v4_basin.bin、diags_set_A.m这类文件真正的主体可能被组织在子目录中。本文按源码包根目录视角来讲解无论原始布局如何都以make.bat所在位置为根目录。2.1 文件清单从 make.bat 到 v4_basin.bin文件作用我的处理方式make.batWindows 下编译 mex 功能的批处理如果后续要用 mex 加速先在 Matlab 中执行mex -setup再双击运行bibli.bib工具箱参考论文的 BibTeX 库写论文时导入即可不参与运行v4_basin.bin盆地掩膜二进制数据用于区分太平洋、印度洋等通常放入数据目录用fread或rdmds读取rdmds.m读取 MITgcm 二进制输出.data/.meta调用入口第 3 章详细讲process2nctiles.m把 face 分块输出整理成 NetCDF tiles用于最终可视化或交给其他语言diags_set_A.m ... diags_set_F.m一组诊断脚本每个脚本生成一类变量先备份再改路径不要被这些文件名误导diags_set_*.m不是需要逐个运行的独立模块而是一组示例脚本每个展示一种诊断模式。我一般会先打开diags_set_A.m看它前 50 行如何把mygrid和状态变量加载进来然后照葫芦画瓢改自己的路径。另一个容易被忽视的是bibli.bib。这份 BibTeX 文件记录的是 gcmfaces 相关方法的原始论文审稿人如果问你的通量计算依据是什么引用这里面的条目远比引用工具箱本身更有说服力。我之前投气候动力学相关期刊时就是靠这份 bib 文件补齐了两个关键引用。2.2 环境配置不要急着点运行% gcmfaces_example_config.m % 在 Octave 上同样适用 gcmfaces_root /path/to/gcmfaces; % 替换为解压根目录 addpath(genpath(gcmfaces_root)); cd(gcmfaces_root); % 检查是否已经有编译好的 mex 文件 if exist(mex/wrappers, dir) addpath([gcmfaces_root /mex]); end这里的关键是把整个目录树都加入 path。rdmds.m内部会调用gcmfaces_global来初始化全局网格对象如果不用genpath运行时会报“未定义函数或变量gcmfaces_global”。cd到根目录的原因在于包里内置了一些相对路径配置文件比如网格信息预设启动时需要从根目录定位这些资源。如果你想把主程序写在工作目录可以在初始化后单独维护一个grid_file_path不必每次都cd回去。常见坑是 Windows 与 Linux 的路径分隔符不一致。Windows 下make.bat能用但如果你直接调用system(make)可能失败因为默认编译器没有配置。Linux/macOS 下make.bat没有用需要手动执行mex -setup再把*.c文件用mex命令逐个编译。使用 Octave 6.4 和 Matlab R2023a 跑同一套脚本语法差异主要在并行部分parfor换成parcellfun即可。2.3 face 数据结构一次切片操作背后的重载当你写出u mygrid.Lon 180时得到的结果不是一个普通逻辑数组而是一个gcmfaces对象内部包含六个面的逻辑矩阵。这就是 gcmfaces 最核心的思路把地球表面切成若干块每个块矩阵尺寸不同但逻辑上是同一个全球场。size(mygrid.Lon)返回一个 cell其中包含每个面的二维尺寸。因此后续所有诊断脚本都不需要你手动拼接 face 边界sum(u.*v)这样的运算会按 face 逐一计算后自动汇总。如果发现计算结果与预期不符先检查mygrid是否被正确初始化——用plot(mygrid.XC, mygrid.YC, .)画一个点云如果看到六个明显分离的扇面并且边界能拼成一个球面说明 face 结构完好。很多人第一次接触 gcmfaces 时总忍不住用squeeze或reshape去处理对象结果破坏了内部元数据导致最终结果全是 NaN这一点要格外注意。3. 用 rdmds.m 与 process2nctiles.m 把 MITgcm 输出变成可用数据MITgcm 的标准输出是二进制.data和附带元信息的.meta二者成对出现。rdmds.m负责读取它们并返回一个gcmfaces变量。但要注意读取后变量内部的 face 维度顺序是face, x, y, z, time而你习惯的矩阵可能是lat, lon, time。所以下一步必须通过process2nctiles.m转成 NetCDF或者自己把 face 提取出来重新排列。3.1 一次完整的读取-拼接流程% 读取第 86400 秒的三维状态场 myState rdmds(state_3d, 86400); % 第二个参数是迭代次数或时间秒数 % 先初始化全局网格如果还没做 gcmfaces_global; % 把 gcmfaces 对象中的第 1 个 face 取出来 % 注意直接索引对象会触发重载推荐这样取 tmp myState.faces{1}; % tmp 是普通 matlab 数组尺寸为 nx x ny x nz size(tmp)这段代码说明一个重要原则不要对gcmfaces对象直接做高维索引例如myState(:,:,1)可能被重载为对所有 face 的同一层做全局切片而不是你想象中的单个 face。先用.faces{1}把原始矩阵取出来再按常规方式处理这也是diags_set_*.m内部常见的写法。rdmds的第二个参数比较特殊在 MITgcm 输出中.meta文件记录了迭代时间或迭代号传入86400时它会自动匹配对应的state_3d.0000008640.data文件匹配规则由文件后缀位数决定。如果文件后缀是 10 位数字要写成0000008640rdmds能自动补零但前提是文件确实存在。% 使用 process2nctiles.m 将 face 拼成 NC 文件的典型前置操作 % 这个脚本内部会把每个 face 单独写出不依赖全局内存 in_file state_3d; out_dir ./nc_out; tile_count 6; z_list 1:50; % 只写前 50 层 % 老版本 gcmfaces 提供 nctiles_write新版本可能用 fwrite2nctiles if exist(fwrite2nctiles, file) 2 fwrite2nctiles(in_file, out_dir, tile_count); else % 常见做法直接用 nctiles_write nctiles_write(in_file, out_dir, z_list); endprocess2nctiles.m不是一个万能转换器它对输入变量名有约定。MITgcm 标准变量Theta、S、U、V可以直接转换如果是用户自定义变量需要先修改脚本里变量名映射表。这也解释了为什么源码包里同时提供diags_set_*.m那组脚本先把诊断结果重命名为标准 MITgcm 变量名再交给process2nctiles转 NC这样下游工具可以直接使用。3.2 参数配置表与检查手段process2nctiles.m开头的几个关键参数如下参数配置值说明file_namestate_3d输入.meta文件的基础名out_dir/data/nc_outNetCDF 输出目录tile_count6face 数量双极网格可能不是 6z_list1:50要转换的深度层减少 I/O 压力overwritetrue是否覆盖已存在的 NC 文件如果你读出来的数据全是 NaN 或面尺寸不匹配优先检查gcmfaces_global中的nFaces和faceDims。MITgcm 输出有 global 和 regional 两种网格gcmfaces 默认预设是 ECCO v4 的 cs32 网格。如果你的模型分辨率不同必须在gcmfaces_global_specs里指定自己的网格文件路径。检查方法很直接% 打印每个面的尺寸与 .meta 文件中记录的 sizes 对比 for f 1:numel(myGrid.faces) fprintf(grid face %d: %d x %d\n, f, size(myGrid.faces{f},1), size(myGrid.faces{f},2)); end for f 1:numel(myState.faces) fprintf(state face %d: %d x %d\n, f, size(myState.faces{f},1), size(myState.faces{f},2)); end如果两个 face 尺寸不一致说明全局网格初始化错了而不是读取错误。此时重新运行gcmfaces_global并确保指定了正确的网格目录。4. diags_set_A/B/C/D/F诊断脚本的配置逻辑与并行批量处理diags_set_*.m这套脚本是整个工具箱里最值得参考的部分也是多数人复制修改的起点。它们不是同一个程序的不同模式而是针对不同输出组合的诊断流程。以下是从文件命名、变量依赖和实际运行推导出的用法以及批量实验时的修改模板。4.1 每个脚本在算什么脚本典型输出依赖输入运行成本diags_set_A.m基本水团温度、盐度、位势密度Theta,S,PHIHYD低diags_set_B.m水平通量U/V 对温度盐度的输运UVEL,VVEL,Theta,S中diags_set_C.m涡度与散度UVEL,VVEL低diags_set_D.m时间平均流场与方差多时间步的 U/V高建议并行diags_set_F.m净表面热通量与淡水通量TFLUX,SFLUX低这些脚本结构基本固定先读取全局网格再用rdmds读取输入字段然后调用 gcmfaces 重载的div、curl、gradient等函数。比如diags_set_C.m里会出现divU div(U, V)。注意这里的div不是 Matlab 自带的divergence它内部处理了 face 边缘的数据交换自己手写最容易出错。4.2 修改路径与变量名一份可直接套用的模板% 修改自 diags_set_A.m 的模板 run_dir /nas/data/RUN_01; % 模型输出目录 grid_dir /nas/data/grid; % 网格文件目录含 v4_basin.bin addpath(genpath(/opt/gcmfaces)); gcmfaces_global; % 加载网格并读取盆地掩膜 mygrid load_grid(grid_dir); % load_grid 的常见实现见下方注释 basin readbin(fullfile(grid_dir, v4_basin.bin), mygrid); % 读取温度单位摄氏度和实际文件一致 Theta rdmds(fullfile(run_dir, state_2d), 86400); Theta convert2gcmfaces(Theta); Theta setCoordSystem(Theta, mygrid); % 用盆地掩膜提取印度洋区域 mask_IO basin 3; % 具体编号需以数据说明为准 Theta_IO Theta .* mask_IO; save(diag_A_IO.mat, Theta_IO, mask_IO);代码中的load_grid和readbin不是 gcmfaces 自带的函数但常见做法是直接在脚本里用freadreshape读取二进制掩膜再通过convert2gcmfaces转成 face 数据。basin 3的编码是我从 ECCO v4 的四个洋盆序号推测的实际编号必须对图确认。重点在于不要相信任何没画过图的盆地号先用pcolor(basin.faces{1})看一眼编号与空间位置是否对应。我曾经因为字节序问题导致印度洋区域出现北极点的温度就是在这里踩的坑。4.3 批量跑多个时间步用 parfor 或 parcellfun如果你有 500 个时间步需要处理逐个运行会非常慢。diags_set_D.m是为批量设计的内部有一个时间循环。这里给出把它改造成并行函数的方法% diag_D_batch.m tic; parfor tid 1:500 if ~exist(sprintf(diag_D_%04d.mat, tid), file) diag_D_single(tid, /nas/data/RUN_01); end end toc; function diag_D_single(tid, run_dir) gcmfaces_global; % 每个 worker 都需要初始化 mygrid getGlobalGrid(); % 自定义函数只读方式获取全局网格 U rdmds(fullfile(run_dir, U), tid); V rdmds(fullfile(run_dir, V), tid); [divUV, curlV] computeDivCurl(U, V, mygrid); save(sprintf(diag_D_%04d.mat, tid), divUV, curlV); endparfor要求每个迭代独立而diags_set_*.m中的全局网格对象如果在新 worker 中没有初始化会报空指针错误。所以我在函数内部显式调用gcmfaces_global用getGlobalGrid从全局变量读取网格避免每次迭代重复加载。Octave 中没有parfor可以用parcellfun但要注意每次迭代结束时输出对象会被序列化到临时文件传输gcmfaces对象时可能会丢失内部缓存。一个稳妥方案是让每个 worker 直接写独立.mat文件最后再统一合并。4.4 常见失败内存溢出与路径分隔符高分辨率 GCM 输出动辄几十 GBrdmds一次性读入全部变量会把内存打爆。我的经验是在运行前先看可用内存是否大于数据文件大小的 1.5 倍如果不够改用diags_set_A.m中注释掉的分块读取模式。另一个高频错误是路径分隔符。diags_set_*.m源文件中如果写死了\到 Linux 下会报“No such file or directory”。把写死的路径改成fullfile组合能同时兼容 Windows 和 Linux。5. 验证 gcmfaces 诊断结果的一个技巧手工温度梯度对照最后一个技巧用不依赖外部工具的方式检验 gcmfaces 算出的“梯度”是否正确。原理很简单构造一个线性温度场用普通有限差分手算梯度再与 gcmfaces 的gradient重载结果对比。如果差异超过数值精度范围说明 face 边界交换或坐标定义出了问题。% verify_gradient.m % 取第 1 个 face 作为规则区域测试 F myGrid.faces{1}; Lon F.XC; Lat F.YC; T 20 0.05 * (Lon - min(Lon(:))) 0.01 * (Lat - min(Lat(:))); % gcmfaces 的梯度重载 dT_gcm gradient(T, myGrid); % 手算经纬向差分 dx mean(diff(Lon(1,:))); % 经向分辨率 dy mean(diff(Lat(:,1))); % 纬向分辨率 [dTdx, dTdy] gradient(T, dx, dy); % 比较差异只比较第 1 个 face 内部 maxDiff max(abs(dT_gcm.faces{1}(:)) - abs(dTdx(:)) ); disp(maxDiff);如果maxDiff接近 0说明 face 内部计算正确。但要特别注意这个测试只能验证单个 face 内部不能验证跨 face 边界。想看跨 face 是否正常需要在两个 face 交界处各取一个点用pcolor画等值线观察是否有错位。我就是在类似测试中抓到过v4_basin.bin读入时字节序错误导致掩膜翻转进而让整个印度洋区域数据异常。这个验证技巧还能用来检查process2nctiles.m输出的 NetCDF 文件。把 NC 文件用ncread读回 Matlab再与 gcmfaces 原生的plot(Theta)对比颜色分布一致才算通过。要注意的是gcmfaces 默认输出网格的坐标顺序不是简单的lon, lat输出 NC 文件时务必把lon、lat数组显式写进变量属性里避免下游 Python or Julia 工具读出来只有 face 编号。最后再补一句任何对盆地掩膜的修改都建议在验证脚本之后再做否则结果做完了才发现掩膜方向反了返工成本极高。本文还有配套的精品资源点击获取