WRF学习笔记之三：使用ERA5数据驱动并运行WRFV4.4（一层嵌套）/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

出的图：

几乎没啥变化（这不是当然吗），毕竟是海温变化不会大。
总之跑通了，至于效果好不好，之后再摸索吧。