1. 从零开始理解MJO位相分析第一次接触MJOMadden-Julian Oscillation马登-朱利安振荡位相分析时我也曾被各种专业术语绕得头晕。简单来说MJO是热带地区一种周期约30-60天的大气波动现象它会显著影响全球天气和气候。而位相分析就是把这种波动按照不同发展阶段分成8个状态就像把一个月亮周期分成新月、上弦月等不同月相。为什么要用NCL和ERA5数据来做这件事NCLNCAR Command Language是气象领域的老牌数据处理工具虽然现在Python越来越流行但NCL在气象专业分析中仍有不可替代的优势。ERA5则是欧洲中期天气预报中心ECMWF提供的第五代再分析数据可以说是目前最全面的全球气候数据集之一。记得我第一次跑这个流程时光是准备数据就花了整整两天。ERA5数据虽然质量高但动辄几十GB的文件下载和处理真是让人头疼。后来发现可以用CDS API批量下载效率提升不少。这里给新手一个建议先从短时间范围比如1年的数据开始练习等流程跑通了再处理更长时间序列。2. 数据准备与环境搭建2.1 获取ERA5数据我们需要四种关键变量OLROutgoing Longwave Radiation向外长波辐射u850850hPa纬向风v850850hPa经向风u200200hPa纬向风在ECMWF官网申请API密钥后可以用Python脚本批量下载。比如下载OLR数据的请求参数长这样{ product_type: reanalysis, variable: toa_outgoing_longwave_radiation, year: 2020, month: [01,02,03], day: [01,02,03], time: [00:00,06:00,12:00,18:00], format: netcdf }2.2 NCL环境配置如果你的系统还没安装NCL推荐用conda安装conda create -n ncl_env -c conda-forge ncl conda activate ncl_env我习惯把脚本和数据按这样的目录结构组织/MJO/ ├── data/ # 存放原始ERA5数据 ├── scripts/ # NCL脚本 ├── output/ # 处理结果 └── plots/ # 生成的图表3. 数据预处理实战3.1 计算异常场原始数据需要先转换成异常场anomalies也就是去掉气候平均态。NCL的clmDayTLL函数可以方便地计算日气候态; 计算日气候态 xClmDay clmDayTLL(x, yyyyddd) ; 用4个谐波平滑 xClmDay_sm smthClmDayTLL(xClmDay, 4) ; 计算异常场 xAnom calcDayAnomTLL(x, yyyyddd, xClmDay_sm)这里有个坑我踩过时间坐标的处理。ERA5数据默认用hours since 1900-01-01这样的格式而NCL脚本可能需要调整时间变量才能正确读取。建议先用ncdump -h查看nc文件的时间维度定义。3.2 数据插值技巧不同变量分辨率可能不一致比如u200是1°×1°而OLR是2.5°×2.5°。我们需要统一插值到相同网格。CDOClimate Data Operators是处理这类任务的神器# 将数据插值到144×181网格 cdo remapbil,r144x181 input.nc output.nc如果遇到Unsupported grid type: generic报错可能是网格定义有问题。可以先提取网格描述文件再修改cdo griddes input.nc mygrid sed -i s/generic/lonlat/g mygrid cdo setgrid,mygrid input.nc output.nc4. EOF分析与MJO信号提取4.1 多变量EOF分析这是整个流程最核心也最难理解的部分。简单说我们要把OLR、u850、u200三个变量的异常场打包在一起做EOF分析找出最主要的空间分布模式。; 合并三个变量的数据 cdata new((3*mlon,ntim), typeof(olr)) do ml0,mlon-1 cdata(ml,:) (/ olr(:,ml) /) cdata(mlmlon,:) (/ u850(:,ml) /) cdata(ml2*mlon,:) (/ u200(:,ml) /) end do ; 计算前两个EOF模态 eof_cdata eofunc(cdata, 2, False)这里有个关键点数据标准化。不同变量的量级差异很大OLR单位是W/m²风场单位是m/s必须先归一化处理olr olr/sqrt(zavg_var_olr) u850 u850/sqrt(zavg_var_u850) u200 u200/sqrt(zavg_var_u200)4.2 位相空间构建得到EOF时间序列后我们可以把每天的MJO状态表示在二维平面上PC1为x轴PC2为y轴。就像钟表的时针位置代表不同时间一样这个平面上的角度位置就对应MJO的不同位相; 计算角度0-360° ang atan2(pc2,pc1)*180./pi nn ind(ang.lt.0) ang(nn) ang(nn) 360 ; 转换为0-360范围 ; 定义8个位相的边界 nPhase 8 angBnd new((/2,nPhase/), float) angBnd(0,:) fspan(0,315,nPhase) ; 0,45,90,...,315 angBnd(1,:) fspan(45,360,nPhase) ; 45,90,...,3605. 位相合成与可视化5.1 位相合成分析最后一步是把每个位相对应的天气场合成起来。这里要注意只选择MJO活跃期PC1²PC2²1; 筛选特定季节和位相的数据 nt ind(mjo_indx.gt.1.0 .and. ; MJO活跃 (imon.ge.5 .and. imon.le.10) .and. ; 5-10月 ang.ge.angBnd(0,n) .and. ang.lt.angBnd(1,n)) ; 位相n ; 计算合成场 xAvg dim_avg_Wrap(x(lat|:,lon|:,time|nt)) uAvg dim_avg_Wrap(u(lat|:,lon|:,time|nt)) vAvg dim_avg_Wrap(v(lat|:,lon|:,time|nt))5.2 专业级可视化NCL的绘图功能非常强大但调图也是门艺术。这是我常用的风场和OLR叠加绘图设置res True resmpFillOn False ; 不填充陆地 rescnFillOn True ; 填充等值线 rescnLinesOn False ; 不画等值线 resgsnScalarContour True ; 矢量标量叠加 resvcRefLengthF 0.04 ; 矢量参考箭头大小 rescnFillPalette ViBlGrWhYeOrRe ; 色标最终生成的位相图应该呈现清晰的东传信号从印度洋位相2-3到海洋大陆位相4-5再到西太平洋位相6-7。如果结果不明显可能需要检查数据时间范围是否足够长至少10年滤波参数是否正确20-100天带通滤波区域平均范围是否合适通常取15°S-15°N6. 常见问题排查在实际操作中我遇到过不少报错和异常情况这里分享几个典型问题的解决方法问题1EOF分析结果异常检查输入数据是否有缺测值确认数据标准化步骤是否正确执行尝试调整经纬度范围热带地区效果最好问题2图形显示异常确保安装了最新版NCL检查色标范围设置是否合理尝试更换图形后端X11转PNG或PDF问题3脚本运行缓慢减少调试时的数据时间范围关闭实时绘图设置PLOTFalse使用NCL的并行计算功能记得第一次成功跑出八位相图时的兴奋感虽然花了整整一周时间调试。现在回头看那些踩过的坑都成了宝贵的经验。建议新手做好两点一是保持耐心气候数据分析本身就是个细致活二是养成写注释的习惯复杂的处理流程隔段时间自己都可能看不懂。