WRF学习笔记之三：使用ERA5数据驱动并运行WRFV4.4（一层嵌套）/WRF运行实录/WRF报错（踩坑）记录
之前写过的和这两篇博文里，我都简单的介绍了WRF的安装与运行，不过，用示例数据来运行只能初步了解WRF的运行流程，想要更进一步熟练地掌握WRF的运行与设置，还需要自己亲自下数据驱动】、设置才行，这篇博文将介绍使用ERA5数据驱动WRF运行，并简单说一下自己犯过的错误，以帮助后学者。注：本次WRF是在组里的小型服务器上的一个结点运行的，所以不涉及作业调度内容。周日晚才拿到的账号，这几天光记着转模式去
·
之前写过的WRF入门和超算入门：大型机运行WRF这两篇博文里，我都简单的介绍了WRF的安装与运行，不过，用示例数据来运行只能初步了解WRF的运行流程，想要更进一步熟练地掌握WRF的运行与设置，还需要自己亲自下数据驱动】、设置才行，这篇博文将介绍使用ERA5数据驱动WRF运行，并简单说一下自己犯过的错误，以帮助后学者。
注：本次WRF是在组里的小型服务器上的一个结点运行的，所以不涉及作业调度内容。
周日晚才拿到的账号，这几天光记着转模式去了，花了一天安装编译，又花了两天才跑通，坑太多了，记录一下。
数据准备
本次使用2018年ERA5数据驱动WRF，下载的是ERA5的single pressure数据和pressure数据，使用python脚本＋IDM批量下载，关于ERA5的下载请参见博文：批量下载ERA5数据（Python+IDM）
使用ERA5驱动数据时有些数据是必须的，因此下载数据这一步骤需要注意不要漏掉变量，
下面是我下载数据的一些代码与信息：
下载气压层（3D数据）
mport cdsapi
import calendar
from subprocess import call
def idmDownloader(task_url, folder_path, file_name):
"""
IDM下载器
:param task_url: 下载任务地址
:param folder_path: 存放文件夹
:param file_name: 文件名
:return:
"""
# IDM安装目录
idm_engine = "C:\\Program Files (x86)\\Internet Download Manager\\IDMan.exe"
# 将任务添加至队列
call([idm_engine, '/d', task_url, '/p', folder_path, '/f', file_name, '/a'])
# 开始任务队列
call([idm_engine, '/s'])
if __name__ == '__main__':
c = cdsapi.Client() # 创建用户
# 数据信息字典
dic = {
'product_type': 'reanalysis', # 产品类型
'format': 'grib', # 数据格式
'variable': [
'geopotential','relative_humidity','specific_humidity',
'temperature','u_component_of_wind','v_component_of_wind'
],
'pressure_level': [
'1','2','3',
'5','7','10',
'20','30','50',
'70', '100','125',
'150', '175',
'200', '225', '250',
'300', '350', '400',
'450', '500', '550',
'600', '650', '700',
'750', '775', '800',
'825', '850', '875',
'900', '925', '950',
'975', '1000',
],
'year': '2018', # 年，设为空
'month': '05', # 月，设为空
'day': [ '01', '02','03','04','05','06','07',], # 日，设为空
'time': [ # 小时
'00:00', '06:00', '12:00',
'18:00'
],
'area': [90, -180, 45, 180],
}
# 通过循环批量下载1979年到2020年所有月份数据
'''
for y in range(1979, 2021): # 遍历年
for m in range(1, 13): # 遍历月
day_num = calendar.monthrange(y, m)[1] # 根据年月，获取当月日数
# 将年、月、日更新至字典中
dic['year'] = str(y)
dic['month'] = str(m).zfill(2)
dic['day'] = [str(d).zfill(2) for d in range(1, day_num + 1)]
'''
r = c.retrieve('reanalysis-era5-pressure-levels', dic, ) # 文件
url = r.location # 获取文件下载地址
path = 'E:\\ERA5\\' # 存放文件夹
filename = 'pressure'+str(2018) + str(5).zfill(2) + '.grib' # 文件名
idmDownloader(url, path, filename) # 添加进IDM中下载
下载单层（2D数据）：
import cdsapi
import calendar
from subprocess import call
def idmDownloader(task_url, folder_path, file_name):
"""
IDM下载器
:param task_url: 下载任务地址
:param folder_path: 存放文件夹
:param file_name: 文件名
:return:
"""
# IDM安装目录
idm_engine = "C:\\Program Files (x86)\\Internet Download Manager\\IDMan.exe"
# 将任务添加至队列
call([idm_engine, '/d', task_url, '/p', folder_path, '/f', file_name, '/a'])
# 开始任务队列
call([idm_engine, '/s'])
if __name__ == '__main__':
c = cdsapi.Client() # 创建用户
# 数据信息字典
dic = {
'product_type': 'reanalysis', # 产品类型
'format': 'grib', # 数据格式
'variable': [
'10m_u_component_of_wind', '10m_v_component_of_wind', '2m_temperature',
'mean_sea_level_pressure', 'sea_surface_temperature', 'skin_temperature',
'snow_albedo', 'snow_depth', 'snowfall',
'soil_temperature_level_1', 'soil_type', 'surface_pressure',
'total_precipitation','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',
],
'year': '2018', # 年，设为空
'month': '05', # 月，设为空
'day': [ '01', '02','03','04','05','06','07',], # 日，设为空
'time': [ # 小时
'00:00', '06:00', '12:00',
'18:00'
],
'area': [90, -180, 45, 180],
}
# 通过循环批量下载1979年到2020年所有月份数据
'''
for y in range(1979, 2021): # 遍历年
for m in range(1, 13): # 遍历月
day_num = calendar.monthrange(y, m)[1] # 根据年月，获取当月日数
# 将年、月、日更新至字典中
dic['year'] = str(y)
dic['month'] = str(m).zfill(2)
dic['day'] = [str(d).zfill(2) for d in range(1, day_num + 1)]
'''
r = c.retrieve('reanalysis-era5-single-levels', dic, ) # 文件
url = r.location # 获取文件下载地址
path = 'E:\\ERA5\\' # 存放文件夹
filename = 'single'+str(2018) + str(5).zfill(2) + '.grib' # 文件名
idmDownloader(url, path, filename) # 添加进IDM中下载
下载后，会出现：
前往相应的路径查找即可。
下好的数据通过FTTP直接传输至服务器，如果需要中转，使用scp命令上传即可。
运行WPS前处理
修改namelist.wps
根据需要修改即可，具体如何修改可参考我的博文：WRF入门，这里贴上我的设置。
&share
wrf_core = 'ARW',
max_dom =1,
start_date = '2018-05-01_12:00:00','2019-09-04_12:00:00',
end_date = '2018-05-07_00:00:00','2019-09-04_12:00:00',
interval_seconds = 21600
/
&geogrid
parent_id = 1,
parent_grid_ratio = 1,
i_parent_start = 1,
j_parent_start = 1,
e_we = 150,
e_sn = 150,
geog_data_res = '10m',
dx = 45000,
dy = 45000,
map_proj = 'polar',
ref_lat = 90.00,
ref_lon = 0.0,
truelat1 = 60.0,
truelat2 = 60.0,
stand_lon = 0.0,
geog_data_path = '/public/home/zhangzl/Build_WRF//WPS_GEOG/'
/
&ungrib
out_format = 'WPS',
prefix = 'FILE',
/
&metgrid
fg_name = 'FILE'
/
~
geogrid
执行/geogrid.exe 如果成功，会出现geo_em.d01文件。
这部分几乎都没有问题，唯一曾出现的问题就是你的静态地理数据下载不全，按照提示一个个下全就好，我个人选择的10m地理数据如下：
其中带有10m和30s的数据一个都不能少。
ungrib
首先创建链接，由于是ERA5数据，Table应为ECMWF：
ln -sf ungrib/Variable_Tables/Vtable.ECMWF Vtable
创建链接后，WPS目录下会出现Table文件。
随后，链接数据：
./link_grib.csh ../ERA5/*.grib
成功后会出现：
两个文件。
随后./ungrib.exe
此处曾出现报错：grib not the same type
原因：输入的文件既有grib1又有grib2 格式，WPS无法同时处理。
解决方法：下载文件时，将文件后缀名写成了.grb，重新下载写成 .grib即可。
运行成功后会出现FILE打头的文件：
metgrid
没有报错，直接运行。/metgrid.exe
成功出现met_em.d开头文件，由于本次只有一层网格，只有d_01。
查看嵌套区域
使用WPS/util/目录下自带的plotgrids.ncl绘制嵌套网格，由于服务器上ncl版本较高，应使用plotgrids_new.ncl绘制。
此时出现报错：
GKS ERROR NUMBER -208 ISSUED FROM SUBROUTINE GOPWK :
--X driver error: DISPLAY environment variable not set
GKS ERROR NUMBER 25 ISSUED FROM SUBROUTINE GESC :
--SPECIFIED WORKSTATION IS NOT OPEN
fatal:Workstation with PID#8 is not open
fatal:Unable to open Workstation-Can't Create
fatal:Unable to access object with id:-4
fatal:PID #-4 can't be found in NhlSetValues
fatal:_NhlCreate:Invalid Parent id #-4
fatal:PID #-4 can't be found in NhlSetValues
fatal:Invalid plot ID=-4 passed to NhlGetBB
fatal:NhlGetValues:PID #-4 is invalid
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 1551 in file $NCARG_ROOT/lib/ncarg/nclscripts/cs m/gsn_code.ncl
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 2704 in file $NCARG_ROOT/lib/ncarg/nclscripts/cs m/gsn_code.ncl
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 10422 in file $NCARG_ROOT/lib/ncarg/nclscripts/c sm/gsn_code.ncl
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 4697 in file $NCARG_ROOT/lib/ncarg/nclscripts/wr f/WRFUserARW.ncl
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 4814 in file $NCARG_ROOT/lib/ncarg/nclscripts/wr f/WRFUserARW.ncl
fatal:["Execute.c":8637]:Execute: Error occurred at or near line 200 in file util/plotgrids_new.ncl
这是由于超算不支持x11格式造成的，打开plotgrids.ncl文件，将type=x11注释，改为type=pdf，即可。
此时WPS目录下出现.pdf文件，使用evince命令查看。由于我在服务器上没有开图形界面的权限，我直接用scp命令把它从服务器上传过来，这不影响。
运行WRF
修改namelist.input
进入WRF/run文件夹，修改namelist.input。
这里报错较多，为了更好的设置相应参数，现在WPS文件夹下：
使用命令：ncdump -h met_em.d01.2018-05-04_12:00:00.nc 查看其中的：
num_st_layers = 4 ;
num_sm_layers = 4 ;
num_metgrid_levels = 38 ;
// global attributes:
:TITLE = "OUTPUT FROM METGRID V4.4" ;
:SIMULATION_START_DATE = "2018-05-04_12:00:00" ;
:WEST-EAST_GRID_DIMENSION = 150 ;
:SOUTH-NORTH_GRID_DIMENSION = 150 ;
:BOTTOM-TOP_GRID_DIMENSION = 38 ;
:WEST-EAST_PATCH_START_UNSTAG = 1 ;
:WEST-EAST_PATCH_END_UNSTAG = 149 ;
:WEST-EAST_PATCH_START_STAG = 1 ;
:WEST-EAST_PATCH_END_STAG = 150 ;
:SOUTH-NORTH_PATCH_START_UNSTAG = 1 ;
:SOUTH-NORTH_PATCH_END_UNSTAG = 149 ;
:SOUTH-NORTH_PATCH_START_STAG = 1 ;
:SOUTH-NORTH_PATCH_END_STAG = 150 ;
:GRIDTYPE = "C" ;
:DX = 45000.f ;
:DY = 45000.f ;
:DYN_OPT = 2 ;
:CEN_LAT = 89.99999f ;
:CEN_LON = 180.f ;
:TRUELAT1 = 60.f ;
:TRUELAT2 = 60.f ;
:MOAD_CEN_LAT = 89.99999f ;
:STAND_LON = 0.f ;
:POLE_LAT = 90.f ;
:POLE_LON = 0.f ;
:corner_lats = 46.77416f, 46.77416f, 46.77416f, 46.77416f, 46.64143f, 46.64143f, 46.64143f, 46.64143f, 46.64143f, 46.64143f, 46.64143f, 46.64143f, 46.50927f, 46.50926f, 46.50926f, 46.50927f ;
:corner_lons = -45.f, -135.f, 135.f, 45.f, -45.1929f, -134.8071f, 134.8071f, 45.19292f, -44.80708f, -135.1929f, 135.1929f, 44.80709f, -45.f, -135.f, 135.f, 45.f ;
:MAP_PROJ = 2 ;
:MMINLU = "MODIFIED_IGBP_MODIS_NOAH" ;
:NUM_LAND_CAT = 21 ;
:ISWATER = 17 ;
:ISLAKE = 21 ;
:ISICE = 15 ;
:ISURBAN = 13 ;
:ISOILWATER = 14 ;
:grid_id = 1 ;
:parent_id = 1 ;
:i_parent_start = 1 ;
:j_parent_start = 1 ;
:i_parent_end = 150 ;
:j_parent_end = 150 ;
:parent_grid_ratio = 1 ;
:sr_x = 1 ;
:sr_y = 1 ;
:NUM_METGRID_SOIL_LEVELS = 4 ;
这部分参数非常重要，在修改namelist.input时，里面的参数需要与读取的met_em.d01.nc文件一致！
time control部分
这部分主要修改起始时间、结束时间与时间步长，根据WPS里的namelist与metgrid生成的met_em,d文件一致即可。
run_days = 5,
run_hours = 0,
run_minutes = 0,
run_seconds = 0,
start_year = 2018,
start_month = 05,
start_day = 01,
start_hour = 12,
end_year = 2018,
end_month = 05,
end_day = 06,
end_hour = 12,
interval_seconds = 21600 %输入文件的间隔时间，单位秒，此处为间隔6小时
input_from_file = .true.,.true., %是否嵌套运行
history_interval = 60, 60, %输出文件间隔，单位为分，此处为每小时输出一个文件
frames_per_outfile = 1, 1, %输出的out文件个数，如果out太大，可以切割输出
restart = .false., %是否是重新运行
restart_interval = 7200, %重新运行out文件的时间间隔，单位分钟
% 以下全为输出文件、输入文件格式，全选2，表示格式为nc
io_form_history = 2
io_form_restart = 2
io_form_input = 2
io_form_boundary = 2
&domains部分
这部分与你WPS的namelist有关，与WPS输出的met_dem文件相关，
time_step = 90,
time_step_fract_num = 0,
time_step_fract_den = 1,
max_dom = 1,
e_we = 150,
e_sn = 150,
e_vert = 40,
dzbot = 10,
dzstretch_s = 1.3,
dzstretch_u = 1.1,
p_top_requested = 5000,
num_metgrid_levels = 38,
num_metgrid_soil_levels = 4,
dx = 30000,
dy = 30000,
grid_id = 1,
parent_id = 1,
i_parent_start = 1,
j_parent_start = 1,
parent_grid_ratio = 1,
parent_time_step_ratio = 1,
feedback = 1,
smooth_option = 0
/
%全部与namelist中保持一致！
physics部分
这部分主要是物理参数化方案的选择，注意查看WRF文档里，参数化方案的一些要求，
&physics
physics_suite = 'CONUS'
mp_physics = 2, %微物理方案，此处选择Lin
cu_physics = 1, %积云参数化方案，此处选择Kain-Fritsch (new Eta) 方案
ra_lw_physics = 1 %长波辐射方案，此处选择rrtn
ra_sw_physics = 5, %短波辐射方案，此处选择Goddard
bl_pbl_physics = 1, %边界层方案，此处选择YSU
sf_sfclay_physics = 1, %地面层方案，应与边界层方案对应，此处为MM5
sf_surface_physics = 2, %地面方案，这里选的Noah
radt = 30, %与你的分辨率对应，dx=20km就写30
bldt = 0, %推荐0
cudt = 0, % 推荐0，除了 Kain-Fritsch参数化方案，这里我忘记改了
icloud = 1,% 是否考虑云辐射效应，这里先不管，直接默认
num_land_cat = 21,%与met_em文件对应
sf_urban_physics = 2,%城市冠层方案，注意与边界层方案对应
fractional_seaice = 1, %把海冰当作分数考虑，如果为0，则要么有冰有么没冰
/
dynamic部分与其他
没什么好讲的……
注意垂直层可能要加个etac上去。
&dynamics
hybrid_opt = 2,
etac = 0.2,
w_damping = 0,
diff_opt = 1,
km_opt = 4,
diff_6th_opt = 0, ,
diff_6th_factor = 0.12, 0.12,
base_temp = 290.
damp_opt = 3,
zdamp = 5000.,
dampcoef = 0.2,
khdif = 0,
kvdif = 0,
non_hydrostatic = .true., .true.,
moist_adv_opt = 1, 1,
scalar_adv_opt = 1, 1,
gwd_opt = 1, 0,
/
&bdy_control
spec_bdy_width = 5,
specified = .true.
/
&grib2
/
&namelist_quilt
nio_tasks_per_group = 0,
nio_groups = 1,
real运行与报错记录
这里报错非常之多，绝大多数都是因为namelist.input设置有误，可以使用：
tail rsl.error.0000 #or tail rsl.out.0000 #or vi rsl.error.0000查看原因并修改。
报错记录一：
You need one of four things:
d01 2018-05-01_12:00:00 1) More eta levels: e_vert
d01 2018-05-01_12:00:00 2) A lower p_top: p_top_requested
d01 2018-05-01_12:00:00 3) Increase the lowest eta thickness: dzbot
d01 2018-05-01_12:00:00 4) Increase the stretching factor: dzstretch_s or dzstretch_u
d01 2018-05-01_12:00:00 All are namelist options
提示需要修改namelist中的四个选项中的一个。
解决方法：查了一下应该输入的数据层数不对应，垂直层间距超过了默认值1000m，因此需要：1、增加垂直层数e_vert;2、降低模式层顶的高度p_top;3、增加允许的最大eta层厚度max_dz。4、增加dz延展因子。
总之就是垂直分辨率不对，按照提示改就行了。
报错记录二：
d01 2018-05-01_12:00:00 t(i,j,k) was 0 at 145 149 30 , setting Qv to 0
d01 2018-05-01_12:00:00 t(i,j,k) was 0 at 146 149 30 , setting Qv to 0
d01 2018-05-01_12:00:00 t(i,j,k) was 0 at 147 149 30 , setting Qv to 0
d01 2018-05-01_12:00:00 t(i,j,k) was 0 at 148 149 30 , setting Qv to 0
d01 2018-05-01_12:00:00 t(i,j,k) was 0 at 149 149 30 , setting Qv to 0
Using sfcprs to compute psfc
-------------- FATAL CALLED ---------------
FATAL CALLED FROM FILE: <stdin> LINE: 6241
troubles, could not find trapping x locations
解决方法：属于网格设置问题，设置的区域超过了数据范围……
重新设置namelist.wps格点，把dx、dy从45km改成30km。
再重新运行WPS。
报错记录三：
The reference pressure is not monotonically decreasing
This tends to be caused by very high topography
(i,j) = 28 11 , topography = 2963.077 m
k = 10 , reference pressure = 67832.36 Pa
k = 11 , reference pressure = 67853.36 Pa
垂直层设置会造成极高的对流层顶，解决方法：不断调整垂直层大小和&dynamic中的etac值（耐心）。
报错记录4：
[zhangzl@node4 run]$ tail rsl.out.0000
column of pressure and value = NaN 12336.15
column of pressure and value = NaN 13345.33
column of pressure and value = NaN 14534.03
column of pressure and value = NaN 15983.12
column of pressure and value = NaN 18277.69
column of pressure and value = -Infinity 20446.78
-------------- FATAL CALLED ---------------
FATAL CALLED FROM FILE: <stdin> LINE: 6241
troubles, could not find trapping x locations
说是垂直气压层上没有数据，很疑惑，明明下载了50hPa的数据，反复修改 层数、etac等依然重复报错。
于是有两个猜想：
1、下载数据时只下到了50hPa，num_metgrid_levels=30太少不够。
2、ungrib链接时将3D数据(pressure)与2D数据(single)同时链接出现问题。
于是重新去下载数据，这次下到了1hPa，将namelist中的 num_metgrid_levels改成38，重新跑ungrib和metgrib，成功运行real，完成了初始化。
wrf.exe运行
继续报错……
报错记录一：
```python
CM_AC_URB3D 0.0000000E+00
SFVENT_URB3D 0.0000000E+00
LFVENT_URB3D 0.0000000E+00
FRC_URB2D 0.0000000E+00
UTYPE_URB2D 0
I 1 J 120
num_urban_hi 15
-------------- FATAL CALLED ---------------
module_physics_init: Use sf_sfclay_physics= 1 or 91 for this pbl option
----------------------------------------
因为我的边界层方案选择YSU，因此 必须sf_sfclay_physics=1。修改继续运行wrf.exe、
可以用ps-a查看进程，一阵漫长的等待后，得wrfout_d01_YYYY-MM_DD_hh:mm:ss 文件。
查看输出结果
用ncl看一下输出文件的变量信息：
Number of global attributes: 134
Number of dimensions: 10
Number of variables: 189
189个变量，有点多，就画下海温吧。
float SST ( Time, south_north, west_east )
FieldType : 104
MemoryOrder : XY
description : SEA SURFACE TEMPERATURE
units : K
stagger :
coordinates : XLONG XLAT XTIM
使用自带NCL自带的脚本简单画一画：
;Load NCL scripts
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
s=addfile("wrfout_d01_2018-05-05_18:00:00","r")
; print(s)
; printVarSummary(s)
sst=wrf_user_getvar(s,"SST_INPUT",0) ;Surface SW downwelling flux
sst@_FillValue =9999
; printVarSummary(dsw)
wks = gsn_open_wks("png","sst_input")
; cmap = read_colormap_file("BlAqGrYeOrRe"); choose colormap
opts = True ; Set some Basic Plot options
opts@MainTitle = "SEA SURFACE TEMPERATURE INPUT" ; plot mods desired
res = opts ; Use basic options for this field
res@cnFillOn = True ; Create a color fill plot
res@ContourParameters = (/275.,280.,1. /) ; Set the levels
contour = wrf_contour(s,wks,sst,res)
pltres = True ; Set plot options
mpres = True ; Set map options
plot = wrf_map_overlays(s,wks,(/contour/),pltres,mpres) ; Plot the field over a map background
end
出的图：
几乎没啥变化（这不是当然吗），毕竟是海温变化不会大。
总之跑通了，至于效果好不好，之后再摸索吧。
华为开发者空间，是为全球开发者打造的专属开发空间，汇聚了华为优质开发资源及工具，致力于让每一位开发者拥有一台云主机，基于华为根生态开发、创新。
更多推荐
16
0- 0
已为社区贡献2条内容
所有评论(0)