Piggy_Packages V2026.1 帮助文档（三）WRF区域模式

获取Piggy_Packages还没有Piggy_Packages的同学请参考这篇帖子获取Piggy_Packages V2026.1 帮助文档一开箱即用在上一节中我们运行了WRF各种理想模式并且体验了气象可视化工具NCVIEW和气象命令行工具CDO。这一节将继续带领你运行WRF区域模式模拟真实的天气。首先我们回到Home目录创建一个名称为Case1的目录并进入。cd mkdir -p Case1 cd Case1新建一个名为FNL的目录并进入。mkdir -p FNL cd FNL下载2026年3月29-31日的FNL气象资料。新建一个脚本文件名叫download_fnl.csh然后打开。touch download_fnl.csh open -e download_fnl.csh将以下内容复制进去然后保存#!/usr/bin/env csh set filelist ( \ https://osdf-director.osg-htc.org/ncar/gdex/d083002/grib2/2026/2026.03/fnl_20260329_00_00.grib2 \ https://osdf-director.osg-htc.org/ncar/gdex/d083002/grib2/2026/2026.03/fnl_20260329_06_00.grib2 \ https://osdf-director.osg-htc.org/ncar/gdex/d083002/grib2/2026/2026.03/fnl_20260329_12_00.grib2 \ https://osdf-director.osg-htc.org/ncar/gdex/d083002/grib2/2026/2026.03/fnl_20260329_18_00.grib2 \ https://osdf-director.osg-htc.org/ncar/gdex/d083002/grib2/2026/2026.03/fnl_20260330_00_00.grib2 \ https://osdf-director.osg-htc.org/ncar/gdex/d083002/grib2/2026/2026.03/fnl_20260330_06_00.grib2 \ https://osdf-director.osg-htc.org/ncar/gdex/d083002/grib2/2026/2026.03/fnl_20260330_12_00.grib2 \ https://osdf-director.osg-htc.org/ncar/gdex/d083002/grib2/2026/2026.03/fnl_20260330_18_00.grib2 \ https://osdf-director.osg-htc.org/ncar/gdex/d083002/grib2/2026/2026.03/fnl_20260331_00_00.grib2 \ ) while($#filelist 0) set syscmd wget --no-check-certificate -N -c $filelist[1] echo $syscmd ... $syscmd shift filelist end运行这个文件下载数据。csh download_fnl.csh* 如遇到网络问题可以参考下方备用方案~~~~~~~~输入以下命令获取备用链接。echo 572R55uY5LiL6L295Zyw5Z2A77yaaHR0cHM6Ly9wYW4uYmFpZHUuY29tL3MvMXFfSXNxNmt5eXpoSVlEWUkzUXpnMWc/cHdkPTEyMzQ | base64 --decode echo 通过备用方式获取资料后打开当前目录open将获取的FNL-Sample.7z拷贝进当前目录然后解压到当前目录。或者在Piggy_Packages中输入以下命令解压7z x FNL-Sample.7z~~~~~~~~用Panoply可以查看下载的FNL文件Panoply fnl_20260331_00_00.grib2~~~~~提示wgrib2是分析处理GRIB格式数据的常用工具例如输入以下命令将数据裁剪到中国经纬度范围wgrib2 fnl_20260331_00_00.grib2 -small_grib 70:140 0:55 output.grib2了解wgrib2更多用法请访问https://www.cpc.ncep.noaa.gov/products/wesley/wgrib2/~~~~~我们先回到上一级目录将WPS模块拷贝到当前目录。cd .. cp /pp_model/WRF-4.7.1/WPS .进入WPS目录修改namelist.wps配置文件。cd WPS open -e namelist.wps清空原始内容将以下内容复制进去sharewrf_core ARW,max_dom 1,start_date 2026-03-29_00:00:00,end_date 2026-03-31_00:00:00,interval_seconds 10800/geogridparent_id 1,parent_grid_ratio 1,i_parent_start 1,j_parent_start 1,e_we 50,e_sn 50,geog_data_res lowres,dx 81000,dy 81000,map_proj lambert,ref_lat 30.0,ref_lon 110.0,truelat1 20.0,truelat2 50.0,stand_lon 108.9,geog_data_path /pp_model/WRF-4.7.1/WPS_GEOG//ungribout_format WPS,prefix FILE,/metgridfg_name FILE/【可选】下载高精度地形数据Piggy_Packages自带低精度地形数据如需使用更高精度的地形数据可以输入如下命令cd /pp_model/WRF-4.7.1/ # 进入WRF目录 ./download_geog.sh # 运行高精度地形数据下载脚本 cd - # 返回刚才的目录下载完成后请将namelist.wps配置文件中的 geog_data_res lowres, 修改为高分辨率选项如 geog_data_res 30s,为了能在较低配置的电脑上运行WRF模式本帮助文档运行的案例分辨率都较低默认无需高精度地形数据。大家可以自己运行更高分辨率的案例。保存后运行绘图命令plotgrids将会生成区域示意图检查无误后我们来准备地形文件。mpirun -np 4 ./geogrid.exegeogrid.exe将会把WPS地形文件插值到我们模拟的区域。运行完成后输入以下命令查看准备的地形文件Panoply geo_em.d01.ncHGT_M是海拔高度变量在运行WRF时应注意避免模拟区域的4个角位于复杂地形上。检查完毕后我们来处理FNL数据。由于FNL数据和GFS同源因此我们使用GFS的Vtable文件来解析FNL数据。ln -s ungrib/Variable_Tables/Vtable.GFS Vtable然后将FNL数据链接到当前目录。./link_grib.csh ../FNL/fnl_*.grib2运行ungrib.exe解码FNL资料然后运行metgrid.exe将FNL资料插值到模式网格上与地形数据合并。./ungrib.exe mpirun -np 4 ./metgrid.exe可以通过NCVIEW或者Panoply查看生成的文件。ncview met_em.d01.2026-03-29_00:00:00.nc返回到上级目录然后将WRF-Real拷贝过来并进入其中的run目录。cd .. cp /pp_model/WRF-4.7.1/WRF-Real . cd WRF-Real/run修改namelist.inputopen -e namelist.input把以下内容贴进去替换原始内容time_controlrun_days 2,run_hours 0,run_minutes 0,run_seconds 0,start_year 2026,start_month 03,start_day 29,start_hour 00,end_year 2026,end_month 03,end_day 31,end_hour 00,interval_seconds 10800input_from_file .true.,history_interval 10,frames_per_outfile 1,restart .false.,restart_interval 7200,io_form_history 2io_form_restart 2io_form_input 2io_form_boundary 2/domainstime_step 120,time_step_fract_num 0,time_step_fract_den 1,max_dom 1,e_we 50,e_sn 50,e_vert 48,dzbot 30.dzstretch_s 1.11dzstretch_u 1.10p_top_requested 5000,num_metgrid_levels 34,num_metgrid_soil_levels 4,dx 81000,dy 81000,grid_id 1,parent_id 0,i_parent_start 1,j_parent_start 1,parent_grid_ratio 1,parent_time_step_ratio 1,feedback 1,smooth_option 0/physicsphysics_suite CONUSmp_physics -1, -1,cu_physics -1, -1,ra_lw_physics -1, -1,ra_sw_physics -1, -1,bl_pbl_physics -1, -1,sf_sfclay_physics -1, -1,sf_surface_physics -1, -1,radt 15, 15,bldt 0, 0,cudt 0, 0,icloud 1,num_land_cat 21,sf_urban_physics 0, 0,fractional_seaice 1,/fdda/dynamicshybrid_opt 2,w_damping 0,diff_opt 2, 2,km_opt 4, 4,diff_6th_opt 0, 0,diff_6th_factor 0.12, 0.12,base_temp 290.damp_opt 3,zdamp 5000., 5000.,dampcoef 0.2, 0.2,khdif 0, 0,kvdif 0, 0,non_hydrostatic .true., .true.,moist_adv_opt 1, 1,scalar_adv_opt 1, 1,gwd_opt 0, 0,/bdy_controlspec_bdy_width 5,specified .true./grib2/namelist_quiltnio_tasks_per_group 0,nio_groups 1,/将上一步metgrid.exe输出的met_em文件链接过来运行real.exe和wrf.exeln -s ../../WPS/met_em.* . mpirun -np 4 ./real.exe mpirun -np 4 ./wrf.exeln -s 是创建软链接的命令一个点号代表当前目录两个点号代表上级目录。这里别忘了第一行命令最后有一个点号代表链接到当前目录。mpirun -np 4 代表用4个进程运行程序。在CPU和内存充足的情况下进程个数越多计算速度通常也会更快但进程个数不可超过CPU核心个数。在wrf运行时可以再打开一个终端查看计算进度。tail -f ~/Case1/WRF-Real/run/rsl.out.0000运行完成后通过Panoply查看wrfout文件Panoply wrfout_d01_2026-03-31_00:00:00下面我们来计算本次模拟过程中2026年3月30日00时至31日00时世界时期间的降水。WRF中常见的降水累计变量有三个RAINC深对流累计降水RAINSH浅对流累计降水RAINNC网格尺度微物理过程产生的累计降水总累计降水通常是RAINC RAINSH RAINNC需要注意的是这里的降水值是累积值即从预报开始到当前时刻的累计降水量。为了获得一段时间的降水量我们需要将两个时刻的累计降水量相减。cdo -expr,precip(RAINCRAINNCRAINSH) -sub wrfout_d01_2026-03-31_00:00:00 wrfout_d01_2026-03-30_00:00:00 precip_24h_20260330.nc命令中-expr 是执行计算的意思-sub是相减的意思。我们用Panoply查看计算结果Panoply precip_24h_20260330.ncWRF输出的垂直层为模式层我们需要通过Python、NCL等工具将它转换为气压层。例如我们来输出500hPa高度的风速新建一个绘图脚本然后打开。touch plot.py open -e plot.py粘贴下方代码import matplotlib.pyplot as plt import numpy as np from netCDF4 import Dataset import cartopy.crs as ccrs import cartopy.feature as cfeature from wrf import getvar, interplevel, latlon_coords, get_cartopy, cartopy_xlim, cartopy_ylim file_path wrfout_d01_2026-03-31_00:00:00 ncfile Dataset(file_path) wspd getvar(ncfile, uvmet_wspd_wdir, unitsm s-1)[0] # 提取风速 p getvar(ncfile, pressure) # 提取气压 wspd_500 interplevel(wspd, p, 500.0) # 插值到500hPa等压面 lats, lons latlon_coords(wspd_500) cart_proj get_cartopy(wspd_500) fig plt.figure(figsize(10, 8)) ax plt.axes(projectioncart_proj) clevs np.arange(0, 60, 5) plt.contourf(lons, lats, wspd_500, levelsclevs, transformccrs.PlateCarree(), cmapjet) plt.colorbar(axax, shrink.8, labelWind Speed (m/s)) ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.LAND, facecolornone, edgecolorblack) ax.add_feature(cfeature.BORDERS, linestyle:) ax.set_xlim(cartopy_xlim(wspd_500)) ax.set_ylim(cartopy_ylim(wspd_500)) gl ax.gridlines(draw_labelsTrue, dmsTrue, x_inlineFalse, y_inlineFalse) gl.top_labels False gl.right_labels False plt.title(500hPa Wind Speed (m/s) - 2026-03-31) plt.savefig(wspd_500hpa_20260331.png, dpi300, bbox_inchestight)打开绘图open wspd_500hpa_20260331.png