用ERA5数据驱动WRFV4.4做区域模拟,这几年几乎成了模式圈里的默认操作。相比GFS和FNL,ERA5的时间跨度更长、资料一致性更好,尤其在上游观测稀疏的地区,做出来的初始场和边界场更干净,跑出来的结果也更稳定。但真正上手时,大家遇到最多的问题往往不是模式本身,而是从数据下载开始就一路踩坑:变量选不对、grib转不出中间文件、real.exe报错、wrf.exe跑几步就崩。
这篇文章直接讲一套能落地的单层嵌套(单域)方案:从CDS下载ERA5、用WPS生成met_em文件、配置WRFV4.4的namelist.input,再到real.exe和wrf.exe的常见报错怎么排查。适合第一次接ERA5跑区域模式的研究生,也适合想从多层嵌套改成单域快速验证方案思路的工程师。我不会把每个参数都念一遍文档,只会讲我实际用过、出过问题并且验证过的配置和排错方式。
1. 项目整体思路与方案选型
1.1 为什么最终选ERA5做输入场
做区域模式的初始场和边界场,常用的无非GFS、FNL、ERA5这几类。GFS和FNL胜在时效高,适合业务预报,但历史数据长度受限于存档策略,而且全球模式在部分区域的偏差会直接影响你跑出来的中小尺度系统。ERA5的优势在于它是通过四维变分同化把大量历史观测融合进模式后得到的再分析产品,水平分辨率约31公里,垂直层从地面到0.1 hPa,时间上能回溯到1950年之后,对研究型任务极其友好。
我在实际项目中最看中的一点是长时间序列的一致性。做气候态对比、个例合成分析时,如果输入场本身在不同年份用了不同版本的全球模式,系统性差异很容易被误认为是模拟的物理信号。ERA5虽然也有版本升级,但整体上更统一,换资料带来的人工跳变要少得多。
另外一个实用原因是,ERA5的GRIB数据可以直接被WPS的ungrib识别,配合WPS自带的Vtable.ERA5,不需要写额外的转换脚本。这比早期处理ECMWF数据时常用的ecmwf_coeffs方案省事不少。后面我会单独讲Vtable的配置,这里是整个流程能不能跑通的第一道关。
1.2 单层嵌套到底解决什么问题
单层嵌套这个概念在WRF里指的就是只跑一个模拟域,namelist里对应max_dom = 1。很多刚接触WRF的人容易陷入“嵌套越多越好”的误区,觉得父域给子域提供边界,子域分辨率更高,结果就一定更细。但嵌套带来的不只是成本翻倍,更重要的是边界处理、反馈机制和横向分辨率比例这些环节都会引入新的误差源。如果没到必须做对流尺度模拟的程度,单域方案完全够用。
比如我做9公里分辨率的区域模拟时,单个域就可以覆盖大部分中尺度系统,不需要父域在几百公里外“供边界”。单层嵌套最大的好处是省时间:下载数据量小、WPS处理快、WRF积分步数少,发现问题也好排查——你在rsl.error里看到的错误基本就是这个域本身的问题,而不是嵌套相互作用导致的。
单域方案的代价是边界条件只能靠ERA5的再分析场不断更新,如果你模拟的区域太小(比如水平格点数不足100×100),边界松弛区的虚假扰动会污染内部结果。这里有个经验值:单域模拟的区域范围最好不要小于1000公里×1000公里,网格距在9到27公里之间比较从容。如果你一定要跑3公里以下的对流尺度,再考虑加嵌套也不迟。
1.3 全流程一览:从CDS下载到wrf.exe落地
整套流程可以拆成四个环节,也是我排查报错时的四个检查点:
- 数据准备:从CDS下载ERA5压力层和单层数据,格式必须选grib
- WPS处理:geogrid生成静态地形数据,ungrib把grib转成中间格式,metgrid把中间格式插值到模式网格
- WRF初始化:real.exe读取met_em文件,生成wrfinput和wrfbdy
- 模式积分:wrf.exe完成实际模拟,输出wrfout
这四个环节中任何一个出问题,后面都跑不动。我排错的时候习惯按这个顺序逐级确认:geo_em是否正常、FILE中间文件是否生成、met_em是否完整、wrfinput是否合理。只要前一级没问题,后一级出错的概率就小很多。下文就按照这个顺序展开。
2. ERA5数据获取与预处理
2.1 压力层和单层数据分别要哪些变量
ERA5数据分两套产品:压力层数据和单层数据。WRF需要高空场和地面场一起才能做完整初始化,所以两套都要下载,缺一不可。
压力层数据我通常下载5个核心变量:位势高度(geopotential)、温度(temperature)、相对湿度(relative_humidity)、U风分量(u_component_of_wind)、V风分量(v_component_of_wind)。垂直层建议覆盖从1000 hPa到100 hPa,一般16层左右就够支撑绝大多数中尺度模拟。如果你研究平流层过程或者想加强模式顶附近的垂直分辨率,可以多下几层到50 hPa甚至10 hPa,但要注意num_metgrid_levels必须和实际下载层数严格对应。
单层数据要复杂一些,至少应包括:表面气压(surface_pressure)、平均海平面气压(mean_sea_level_pressure)、2米温度(2m_temperature)、2米露点温度(2m_dewpoint_temperature)、10米U风分量、10米V风分量、地表温度(skin_temperature),以及土壤温度和土壤湿度各四层。如果研究海气过程,海表温度也要下。雪相关的变量在冬季或有积雪覆盖的区域需要注意,建议把雪深、雪水当量、雪密度一起选上。
两个下载时最容易出的问题:一是只下了压力层忘了单层,结果metgrid后met_em里缺TSK、PSFC等关键变量,real.exe直接报错;二是在CDS里选了NetCDF格式而不是grib,导致ungrib认不出文件。我的建议很简单,一律用grib格式。
2.2 cdsapi下载脚本与参数细节
用Python的cdsapi包下载是最方便的方式。先安装并配置CDS API密钥:在用户根目录下创建~/.cdsapirc文件,写入url和key,或者用环境变量配置。下载脚本的核心逻辑是同时请求压力和单层两个数据集,然后按模拟时间段循环。
下面这个脚本是我常用的模板,覆盖了2天时长的9公里分辨率单域模拟所需数据:
import cdsapi c = cdsapi.Client() years = ['2022'] months = ['07'] days = ['01', '02'] times = ['00:00', '06:00', '12:00', '18:00'] area = [45, 100, 25, 125] # 北、西、南、东,按目标区域裁剪 # 压力层数据 c.retrieve( 'reanalysis-era5-pressure-levels', { 'product_type': 'reanalysis', 'variable': [ 'geopotential', 'relative_humidity', 'temperature', 'u_component_of_wind', 'v_component_of_wind' ], 'pressure_level': [ '1000', '975', '950', '925', '900', '850', '800', '700', '600', '500', '400', '300', '250', '200', '150', '100' ], 'year': years, 'month': months, 'day': days, 'time': times, 'area': area, 'data_format': 'grib' }, 'era5_pl_202207.grib' ) # 单层数据 c.retrieve( 'reanalysis-era5-single-levels', { 'product_type': 'reanalysis', 'variable': [ '10m_u_component_of_wind', '10m_v_component_of_wind', '2m_dewpoint_temperature', '2m_temperature', 'mean_sea_level_pressure', 'sea_surface_temperature', 'skin_temperature', 'surface_pressure', 'soil_temperature_level_1', 'soil_temperature_level_2', 'soil_temperature_level_3', 'soil_temperature_level_4', 'volumetric_soil_water_layer_1', 'volumetric_soil_water_layer_2', 'volumetric_soil_water_layer_3', 'volumetric_soil_water_layer_4', 'snow_albedo', 'snow_density', 'snow_depth', 'snow_depth_water_equivalent', 'snowfall' ], 'year': years, 'month': months, 'day': days, 'time': times, 'area': area, 'data_format': 'grib' }, 'era5_sl_202207.grib' )这里area参数的顺序一定要按北、西、南、东来,写反了会下载到完全不对的区域。另外,CDS在高峰期经常排队,实测下载2天、16层、覆盖10°×20°区域的数据量大约在4到8GB左右,视变量数量决定。如果网络不稳定,可以把请求拆小,按天或者按变量组下载,再手动整理到同一个目录。
2.3 区域裁剪与时间拼接
CDS的area参数其实已经做了区域裁剪,这一步可以不用额外处理。但有个细节值得注意:如果你的模拟域中心在东西经边界附近,比如处理跨太平洋区域,就要注意ERA5的经度范围是0到360度,不能用习惯的-180到180度直接填。我遇到过同事用area=[40, -120, 20, -100]下载,结果返回数据为空的案例,后来把西经换算成东经(-120度对应240度)才正常。
如果模拟时间较长,建议下载时按年或按月设置循环,避免单次请求数据量过大导致CDS任务中途失败。下载完成后,可以用grib_ls快速检查文件里的记录数和变量数:
grib_ls -p parameterName,shortName,level,dataDate,dataTime era5_pl_202207.grib这个命令能帮你确认每个时次都有完整的变量和层次。如果发现某个时次缺变量,不要心存侥幸,直接重新下载缺失的时间段,否则ungrib阶段会出现奇怪的缺失字段错误。
2.4 下载和处理过程中的典型坑
第一个坑是变量名写错。CDS的变量名是固定的,比如“10米U风分量”必须写成10m_u_component_of_wind,少个m都不行。建议从CDS官网的变量列表里直接复制,不要手打。
第二个坑是时间间隔不一致。下载间隔决定了后续interval_seconds的设置。我一般下载6小时间隔的数据,对应WRF里interval_seconds = 21600。如果你下载了3小时或1小时间隔的数据,interval_seconds也要同步修改,否则real.exe在生成侧边界时会找不到匹配时次。
第三个坑是GRIB格式和NetCDF格式的选择。WPS的ungrib只能直接处理grib格式,如果选了NetCDF,还要额外写脚本转成WPS中间格式,完全没有必要。直接在data_format参数里写grib就行。
3. WPS:从grib到met_em的关键三步
3.1 geogrid:静态地理数据检查
WPS的第一步是geogrid,它根据你在namelist.wps里定义的模拟域生成geo_em文件,包含地形高度、土地利用类型、土壤类型等静态数据。这些数据不来自ERA5,而是来自GEOG静态数据集。
geogrid配置的关键是参数一致性。如果后面要跑单域,namelist.wps里只需要定义一套域参数:
&geogrid parent_grid_ratio = 1, i_parent_start = 1, j_parent_start = 1, e_we = 160, e_sn = 120, geog_data_res = 'default', dx = 9000, dy = 9000, map_proj = 'lambert', ref_lat = 34.0, ref_lon = 108.0, truelat1 = 30.0, truelat2 = 60.0, stand_lon = 108.0, geog_data_path = '/path/to/WPS_GEOG' /这里e_we和e_sn是网格格点数,我一般设置160×120,对应约1440公里×1080公里的范围,配合9公里分辨率比较合适。geog_data_path必须指向解压后的WPS_GEOG目录,而且目录里要能看到geogrid子目录,否则会报找不到静态数据类型。
geogrid跑完后,用ncdump -h geo_em.d01.nc看一下变量的维度和范围,确认域坐标不是反的、地形高度量级正常。常见问题是经纬度范围定义错误导致模拟域跑到海上或完全不在地球上。检查确认无误后再进入ungrib。
3.2 ungrib:链接顺序与Vtable.ERA5
ungrib的作用是把GRIB格式的原始气象数据解压成WPS中间格式(前缀为FILE的文件)。这一步最关键的坑在Vtable。WPS编译完成后,默认的Vtable链接可能指向通用数据表,而ERA5必须用专门的Vtable.ERA5。
我的做法是进入WPS主目录后,先删掉旧的Vtable链接,再重新建立:
cd WPS rm -f Vtable ln -sf ungrib/Variable_Tables/Vtable.ERA5 Vtable接下来用link_grib脚本把你的grib文件做成ungrib能识别的序列。ungrib通过文件名的后缀顺序读取,所以压力层和单层数据要一起链接进来:
./link_grib.csh /path/to/era5_pl_202207.grib /path/to/era5_sl_202207.grib执行后会生成一堆GRIBFILE.AAA、GRIBFILE.AAB这样的链接。这里有个经验排序:如果你有多个时间段多个文件,建议按时间顺序链接,确保同一天的pressure levels和single levels相邻,这样ungrib能把高空气象和地面气象合并到同一个中间文件里。
然后运行ungrib:
./ungrib.exe正常情况会生成FILE:2022-07-01_00这样的中间文件,文件名里的日期来自grib数据本身。如果生成不了文件,最可能的原因就是Vtable没链接对,或者grib文件的变量与Vtable不匹配。此时回看ungrib的输出,会有一长串“找不到参数”之类的提示。
3.3 metgrid:中间文件的生成验证
metgrid负责把ungrib生成的FILE中间文件插值到geo_em定义的网格上,输出met_em文件。namelist.wps里相应的配置如下:
&metgrid fg_name = 'FILE', io_form_metgrid = 2 /运行前必须确保namelist.wps中的模拟域参数与geogrid完全一致,尤其e_we、e_sn、dx、dy四个值,任何一处不一致都会导致metgrid报“domain size mismatch”之类的错误。
运行:
./metgrid.exe结束后应当看到met_em.d01.2022-07-01_00:00:00.nc这一类文件。我建议用ncdump检查一下met_em里的TSK、PSFC、HGT等变量,确认没有缺变量或量级异常。有一回我下载单层数据时漏了表面气压,met_em里的PSFC全是0,run real.exe时直接报错提示“missing pressure”,排查半天才发现是数据源的问题。
4. 单层嵌套的WRF核心配置
4.1 namelist.input中最容易被忽略的参数
WRF的配置文件是namelist.input,单域模式下的核心参数比你想的简单,但有几个地方特别容易出错。
第一是time_step。它和水平分辨率直接相关,一般经验公式是time_step ≤ 6 × dx(km),也就是9公里分辨率建议用54秒或60秒。用90秒虽然也能跑,但在复杂地形区域容易出现CFL条件超限。我第一次跑9公里域用了90秒,结果在山区附近不断报“cfl”警告,后来降到60秒才稳定。
第二是e_vert和p_top_requested组合。垂直层数要和num_metgrid_levels匹配,而模式层顶不要超过ERA5数据提供的压力层顶部太多。我下载到100 hPa,模式层顶一般设在50 hPa也就是5000 Pa,这样高层有足够的余量让模式自己光滑过渡。
第三是interval_seconds,前面说过,必须和下载的时间间隔一致。
单域namelist.input的核心部分可以这样写:
&time_control run_days = 2, run_hours = 0, start_year = 2022, start_month = 07, start_day = 01, start_hour = 00, end_year = 2022, end_month = 07, end_day = 03, end_hour = 00, interval_seconds = 21600, input_from_file = .true., history_interval = 60, frames_per_outfile = 1, restart = .false., io_form_history = 2, io_form_restart = 2, io_form_input = 2, io_form_boundary = 2, debug_level = 0 /&domains time_step = 60, max_dom = 1, e_we = 160, e_sn = 120, e_vert = 41, p_top_requested = 5000, num_metgrid_levels = 16, dx = 9000, dy = 9000, numtiles = 8 /num_metgrid_levels设为16,这个数字是你在CDS里下的压力层层数,不能随意填。
4.2 物理方案选择:分辨率对应的搭配逻辑
物理方案的选择直接影响模拟结果,而且和分辨率强相关。9公里分辨率属于“灰色地带”,既不是完全分辨对流,也不能简单忽略对流参数化。我实测下来的一套组合如下:
- 微物理:
mp_physics = 8(Thompson方案),对中尺度降水表现稳定 - 长波辐射:
ra_lw_physics = 4(RRTMG) - 短波辐射:
ra_sw_physics = 4(RRTMG) - 近地面层:
sf_sfclay_physics = 2 - 陆面过程:
sf_surface_physics = 2(Noah),单域跑2天足够 - 边界层:
bl_pbl_physics = 6(YSU) - 积云参数化:
cu_physics = 3(Grell-Freitas)
对应的物理配置段:
&physics mp_physics = 8, ra_lw_physics = 4, ra_sw_physics = 4, radt = 15, sf_sfclay_physics = 2, sf_surface_physics = 2, bl_pbl_physics = 6, cu_physics = 3, cu_rad_feedback = .true. /这套组合在夏季降水个例里的表现比较稳。如果你模拟的是冬季强冷空气过程,可以考虑把微物理换成mp_physics = 10(Morrison)或者mp_physics = 9(Milbrandt-Yau),看你对雪和冰相过程的敏感度。但第一次跑通流程时建议先用默认组合,跑通后再换方案对比,不要一上来就做物理敏感性实验。
4.3 垂直层与积分时长设置
垂直层数不需要追求多。41层配合16层的ERA5输入场,已经能较好地刻画边界层和中层大气结构。盲目加到60层或者80层,一方面会增加计算量,另一方面对初始场的垂直插值要求也更高,反而会增加顶层振荡的风险。
积分时长方面,单域跑2到3天是比较合理的验证周期。如果跑7天以上,边界场每6小时更新一次,长时间积分下边界松驰区的误差会逐渐向域内传播,需要额外关注spec_bdy_width等边界设置。单域情况下侧边界条件设置为:
&bdy_control spec_bdy_width = 5, spec_zone = 1, relax_zone = 4, specified = .true., nested = .false. /这里的specified = .true.表示使用指定边界,单域模拟必须这么设。如果你不小心写成了.false.,real.exe不会报错,但wrf.exe运行时边界会变成自由边界,模拟结果会迅速失真。
5. 运行报错排查:从real.exe到wrf.exe
5.1 mpi启动方式与运行队列
WPS的三个程序(geogrid、ungrib、metgrid)通常单进程运行就够了,但real.exe和wrf.exe需要并行。我用的是Intel MPI,启动命令不复杂:
mpirun -np 32 ./real.exe mpirun -np 32 ./wrf.exe关键是在运行real.exe之前,先确认namelist.input里的所有参数与namelist.wps一致,并且met_em文件的时次覆盖了你设定的起止时间。如果met_em文件缺失某个时次,real.exe会卡在找不到输入文件的报错上。
跑完real.exe后,立刻检查是否生成了wrfinput_d01和wrfbdy_d01两个文件。如果只有wrfinput没有wrfbdy,说明侧边界生成失败,这时可以先看rsl.out.0000和rsl.error.0000,一般在最后几行就能看到具体原因。
wrfbdy文件的大小也能反映问题。正常2天模拟、6小时更新一次边界,wrfbdy应该在几百MB量级;如果只有几十KB,那基本等于边界没生成好,wrf.exe跑起来大概率会崩。
5.2 real.exe常见错误定位
real.exe的错误集中在输入场和模式配置不匹配上。我遇到最多的是下面几类。
第一类是关于“找不到met_em文件”。多半是metgrid没跑或者输出路径不对。确认met_em.d01.*文件在WPS目录下或者你指定的目录下存在,并且文件名里的日期和时间段覆盖模拟起止时间。
第二类是“num_metgrid_levels与实际数据不一致”。如果你下载的是16层压力数据,但namelist.input里写成了37,real.exe会在读取垂直剖面时报范围越界。这里没有玄学,下载了几层就写几层。
第三类是地形差异过大。real.exe在插值初始场时会把met_em里的地形高度和geo_em里的地形高度做比较,如果差异超过一定阈值,会输出“inconsistency”类的警告。ERA5的地形分辨率较低,和9公里网格地理数据有一定差异是正常的,但如果差异超过几百米,就要检查是不是模拟域坐标系设置错了。
5.3 wrf.exe常见错误定位
wrf.exe的错误最让人头疼,因为它有时候不是立刻挂掉,而是跑了几百步之后突然出现浮点溢出。排查这类问题,我先看两个文件:rsl.out.0000和rsl.error.0000。
最常见的崩溃原因是CFL条件超限。表现是在rsl.error里反复出现“cfl”字样,然后某个进程报浮点异常退出。这类问题通常和time_step过大、地形过于陡峭,或者物理方案在极端气象条件下不稳定有关。我的处理顺序:先把time_step从60秒降到45秒再试;如果还崩,检查域内是否包含特别高的山脉,考虑把p_top_requested降低;最后再考虑换物理方案。
第二个常见问题是积分过程中出现NaN。这类问题要看NaN最早出现在哪个变量上。如果出现在湿度变量,多半是微物理方案的问题,可以考虑换mp_physics;如果出现在温度变量,要检查辐射方案是否和相关配置冲突。有一个官方文档里不起眼但很关键的参数是damp_opt和zdamp,高层阻尼没开好的时候,模式层顶附近的波动会往上堆积,最终导致整个模式崩溃。
第三个常见问题是并行效率异常或者MPI进程直接中断。这种问题优先检查环境变量,比如ulimit -s unlimited有没有设置,Intel MPI和OpenMPI的混用会不会造成冲突。我自己遇到过OpenMPI编译的WRF用Intel MPI启动后,每个进程都报“unrecognized option”的情况,最后统一了编译器才解决。
5.4 报错速查表
下面这个表是我实际排错过程中总结的高频问题,基本覆盖了从ERA5下载到wrf.exe崩溃的绝大多数场景:
| 现象 | 可能原因 | 解决办法 |
|---|---|---|
| ungrib提示找不到变量或字段 | Vtable没指向Vtable.ERA5,或grib文件下载不完整 | 检查Vtable链接,重新下载缺失数据 |
| ungrib报“corrupted attribute list” | 输入不是grib格式,或者文件被截断 | 确认data_format='grib',重新下载文件 |
| metgrid生成不了met_em文件 | 中间文件前缀不对,或者fg_name与FILE前缀不一致 | 检查ls FILE*,确认fg_name='FILE' |
| real.exe报找不到met_em | 时间范围对不上,或metgrid未成功输出 | 检查起止时间和met_em文件列表 |
| real.exe报垂直层数错误 | num_metgrid_levels和实际下载层数不一致 | 修改namelist.input中的num_metgrid_levels |
| wrf.exe跑几步就cfl | time_step过大 | 将time_step降到45秒或60秒 |
| wrf输出NaN | 物理方案组合不合适 | 检查哪个变量先出现NaN,更换对应方案 |
| MPI进程总是不正常退出 | 环境变量或编译器不匹配 | 统一编译器和MPI类型,设置ulimit -s unlimited |
| wrfout里高层风场异常大 | 模式层顶附近阻尼不足 | 调整damp_opt、zdamp、dampcoef |
这张表不能覆盖所有可能,但能帮你快速定位80%的问题。剩下的20%,基本都能在rsl.out.0000里的警告和提示中找到线索。
6. 一些踩过坑后的经验总结
说实话,用ERA5跑WRF最花时间的不是WRF本身,而是数据准备和排查那些看起来莫名其妙的小错。我自己的习惯是每次改动配置前先备份一份能正常跑通的namelist,同时在WPS目录里保留一份不覆盖的namelist.wps.bak,这样一旦新配置跑挂,可以快速回退对比,不用从头猜。
还有一个容易被忽略的点:诊断输出一定要开够。namelist.input里的history_interval可以设小一点,比如60分钟甚至30分钟,虽然占磁盘,但出问题的时候你至少能知道是哪个时刻、哪个物理量开始异常。很多wrf.exe的崩溃不是一开始就崩,而是某个物理量在边界区累积误差后才崩,没有足够密的输出很难定位。
最后分享一个小技巧:在运行wrf.exe之前,可以用ncdump -h wrfinput_d01快速检查一下初始场的基本量级。如果PSFC的量级在1000 hPa附近、HGT地形高度和你的geo_em一致,那基本可以放心跑。如果初始场就有问题,跑再久也是错的,不要浪费时间。