高性能计算服务
>
软件专区
>
最佳实践
>
WRF台风模拟地形敏感性实验教程，以“格美”为例
WRF台风模拟地形敏感性实验教程，以“格美”为例
WRF(Weather Research and Forecasting)数值模式是当前最广泛使用的中尺度气象模拟系统，其通过嵌套网格技术显著提升低空气象要素的模拟精度。
本次实操，我们以 2024 年台风“格美”为例，基于WRF（Weather Research and Forecasting）模式，结合 NCL 工具，系统介绍如何设计并执行一套完整的地形敏感性实验。实验通过修改地形数据来研究地形对台风路径影响。
实验通过地形敏感性分析，验证了WRF模式在模拟台风与地形相互作用方面的能力，并强调了地形数据精度在台风预报中的重要性。
一、实验简介
本教程将指导您在超算互联网平台使用WRF模式对台风"格美"进行地形敏感性实验，通过修改地形数据来研究地形对台风路径影响。实验将使用NCL（NCAR Command Language）工具修改地形数据，并比较不同地形条件下的台风模拟结果。
二、实验准备
1、环境配置
首先在超算互联网平台搜索wrf4.5,点击购买，点击去使用。 
通过命令行快速打开软件。 
首先加载必要的环境模块并设置Python环境：
shell
module use /public/software/modules/base /public/software/modules/apps
source /work/home/username/apprepo/ncl/6.6.2-gcc4.8.5/scripts/env.sh
source /work/home/username/apprepo/anaconda3/2024.02.1-none/scripts/env.sh
conda create -n wrfdeal python==3.9
conda activate wrfdeal
conda install -c conda-forge wrf-python
pip install netcdf4
pip install meteva2、数据准备
下载FNL初始场数据
使用以下Python脚本下载FNL数据：
shell
import sys
import os
from datetime import datetime, timedelta
from urllib.request import urlretrieve
def generate_gdas1_urls(start_date, end_date):
base_url = "https://data.rda.ucar.edu/d083002/grib2/{year}/{year}.{month:02d}/fnl_{year}{month:02d}{day:02d}_{hour:02d}_00.grib2"
urls = []
current_date = start_date
while current_date <= end_date:
if current_date.hour in [0, 6, 12, 18]:
url = base_url.format(
year=current_date.year,
month=current_date.month,
day=current_date.day,
hour=current_date.hour
)
urls.append(url)
print(url)
current_date += timedelta(hours=6)
return urls
def download_files(urls):
for url in urls:
filename = os.path.basename(url)
print(f"Downloading {filename} ... ", end='', flush=True)
try:
urlretrieve(url, filename)
print("done")
except Exception as e:
print(f"Failed: {e}")
if __name__ == "__main__":
if len(sys.argv) != 3:
print("Usage: python download_gdas1.py <start_date> <end_date>")
print("Date format should be YYYYMMDDHH, for example: 2008072512")
sys.exit(1)
start_str = sys.argv[1]
end_str = sys.argv[2]
try:
start_date = datetime.strptime(start_str, "%Y%m%d%H")
end_date = datetime.strptime(end_str, "%Y%m%d%H")
if start_date > end_date:
print("Error: Start date must not be later than end date.")
sys.exit(1)
file_urls = generate_gdas1_urls(start_date, end_date)
download_files(file_urls)
except ValueError as e:
print(f"Invalid date format: {e}")
print("Date format should be YYYYMMDDHH, for example: 2008072512")
sys.exit(1)运行脚本下载数据：
shell
python download_gdas1.py 2008072512 2008073012三、实验流程
1. 创建项目目录
shell
mkdir -p Project/WRF Project/WPS
cd Project2. 准备WRF和WPS
shell
cp -rf /work/home/username/apprepo/wrf_wps/4.5-intelmpi2017/app/WPS/ WRF
cp -rf /work/home/username/apprepo/wrf_wps/4.5-intelmpi2017/app/WRF/ WPS3. 修改地形数据
使用脚本修改地形数据：
shell
import netCDF4 as nc
import numpy as np
import wrf
# 1. 打开 geo_dem.d01
file_dir = '/work/home/acin2esoa8/project/WPS/geo_em.d01.nc'
data = nc.Dataset(file_dir, 'r+')
# 2. 读取经纬度变量（XLAT_M 和 XLONG_M）
lat = data['XLAT_M'][0, :, :] # 纬度（2D数组）
lon = data['XLONG_M'][0, :, :] # 经度（2D数组）
# 3. 定义台湾省的经纬度范围（大致范围）
taiwan_lat_min, taiwan_lat_max = 21.5, 25.5 # 台湾纬度范围
taiwan_lon_min, taiwan_lon_max = 119.0, 123.0 # 台湾经度范围
# 4. 找到台湾省范围内的网格点（生成一个布尔掩码）
mask = (
(lat >= taiwan_lat_min) & (lat <= taiwan_lat_max) &
(lon >= taiwan_lon_min) & (lon <= taiwan_lon_max)
)
# 5. 修改 HGT_M 的值（示例：台湾省地形高度降低 30%）
hgt = data['HGT_M'][:] # 读取原始地形数据
hgt[0, mask] *= 0.5 # 仅修改台湾省范围内的地形
data['HGT_M'][:] = hgt # 写回文件
# 6. 关闭文件
data.close()
print("台湾省地形高度修改完成！")运行python脚本：
shell
Python ter_mod.py地形修改对比如下： 
4. 运行WPS
修改namelist.wps
shell
cd WPS
vi namelist.wps示例配置：
shell
&share
wrf_core = 'ARW',
max_dom = 1,
start_date = '2024-07-22_00:00:00',
end_date = '2024-07-25_18:00:00',
interval_seconds = 10800
/
&geogrid
parent_id = 1,
parent_grid_ratio = 1,
i_parent_start = 1,
j_parent_start = 1,
e_we = 260,
e_sn = 211,
geog_data_res = 'default',
dx = 9000,
dy = 9000,
map_proj = 'lambert',
ref_lat = 21.00,
ref_lon = 128.00,
truelat1 = 20.0,
truelat2 = 20.0,
stand_lon = 129.6,
geog_data_path = '/public/software/meteorology/WPS_GEOG'
/
&ungrib
out_format = 'WPS',
prefix = 'FILE',
/
&metgrid
fg_name = 'FILE'
/运行WPS预处理：
shell
./geogrid.exe
ln -sf ungrib/Variable_Tables/Vtable.GFS Vtable
./link_grib.csh ../fnl_*
./ungrib.exe
./metgrid.exe5. 运行WRF
修改namelist.input
shell
cd ../WRF/run
vi namelist.input示例配置：
shell
&time_control
start_year = 2024, 2019,
start_month = 07, 09,
start_day = 22, 04,
start_hour = 00, 12,
end_year = 2024, 2019,
end_month = 07, 09,
end_day = 25, 06,
end_hour = 18, 00,
interval_seconds = 10800
input_from_file = .true.,.true.,
history_interval = 60, 60,
frames_per_outfile = 1, 1,
restart = .false.,
restart_interval = 7200,
io_form_history = 2
io_form_restart = 2
io_form_input = 2
io_form_boundary = 2
/
&domains
time_step = 15,
time_step_fract_num = 0,
time_step_fract_den = 1,
max_dom = 1,
e_we = 260, 220,
e_sn = 211, 214,
e_vert = 45, 45,
dzstretch_s = 1.1
p_top_requested = 5000,
num_metgrid_levels = 34,
num_metgrid_soil_levels = 4,
dx = 9000,
dy = 9000,
grid_id = 1, 2,
parent_id = 0, 1,
i_parent_start = 1, 53,
j_parent_start = 1, 25,
parent_grid_ratio = 1, 3,
parent_time_step_ratio = 1, 3,
feedback = 1,
smooth_option = 0
/
&physics
mp_physics = 6, -1,
cu_physics = 0, -1,
ra_lw_physics = 1, -1,
ra_sw_physics = 1, -1,
bl_pbl_physics = 1, -1,
sf_sfclay_physics = 1, -1,
sf_surface_physics = 2, -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
/
&dynamics
hybrid_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 = 1, 0,
/
&bdy_control
spec_bdy_width = 5,
specified = .true.
/
&grib2
/
&namelist_quilt
nio_tasks_per_group = 0,
nio_groups = 1,
/运行WRF：
shell
ln -sf ../../WPS/met_em.d* .
./real.exe
# 提交作业
cp /work/home/username/apprepo/wrf_wps/4.5-intelmpi2017/case/wrf.slurm .
sbatch wrf.slurm6.控制实验
重复上述步骤，不修改地形即可
四、结果分析
shell
使用Python脚本分析不同地形条件下的台风路径和强度：
import numpy as np
from netCDF4 import Dataset
from wrf import (to_np, getvar, smooth2d, get_cartopy, cartopy_xlim,
cartopy_ylim, latlon_coords)
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import glob
# 定义绘图函数
def plot_typhoon_track(wrf_files, titles):
plt.figure(figsize=(15, 10))
ax = plt.axes(projection=ccrs.PlateCarree())
colors = ['red', 'green'] # Only two colors now
for i, (files, title) in enumerate(zip(wrf_files, titles)):
lons, lats, slps = [], [], []
for file in sorted(files):
ncfile = Dataset(file)
slp = getvar(ncfile, "slp")
smooth_slp = smooth2d(slp, 3, cenweight=4)
# 找到最低气压点
min_idx = np.unravel_index(np.argmin(to_np(smooth_slp)), smooth_slp.shape)
lon = to_np(ncfile['XLONG'][0,:,:])[min_idx]
lat = to_np(ncfile['XLAT'][0,:,:])[min_idx]
lons.append(lon)
lats.append(lat)
slps.append(to_np(smooth_slp)[min_idx])
# 绘制路径
ax.plot(lons, lats, color=colors[i], marker='o', linestyle='-',
linewidth=2, markersize=6, label=title, transform=ccrs.PlateCarree())
# 添加地图要素
ax.coastlines(resolution='10m')
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
ax.set_extent([115, 130, 15, 30])
plt.title('Typhoon Track Comparison with Different Terrain Conditions')
plt.legend()
plt.savefig('terrain_sensitivity_track.png', dpi=300, bbox_inches='tight')
plt.close()
# 使用glob自动查找文件
orig_files = sorted(glob.glob('wrfout_orig_d01_*')) # 原始地形结果文件
half_files = sorted(glob.glob('wrfout_half_d01_*')) # 减半地形结果文件
# 绘图比较 (only two scenarios now)
plot_typhoon_track([orig_files, half_files],
['Original Terrain', 'Half-height Terrain'])
import numpy as np
from netCDF4 import Dataset
from wrf import (to_np, getvar, smooth2d)
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
from matplotlib.colors import ListedColormap
from meteva.base.tool.plot_tools import add_china_map_2basemap
def setup_map(ax, extent):
"""Setup map with common elements using meteva's China map"""
# Add China map layers
add_china_map_2basemap(ax, name="river", edgecolor='k', lw=0.5, encoding='gbk') # Rivers
add_china_map_2basemap(ax, name="nation", edgecolor='k', lw=0.5, encoding='gbk') # National borders
add_china_map_2basemap(ax, name="province", edgecolor='k', lw=0.5, encoding='gbk') # Provincial borders
# Set extent
ax.set_extent(extent, crs=ccrs.PlateCarree())
# Configure gridlines
gl = ax.gridlines(draw_labels=True, linewidth=1, color='none', alpha=0.5,
linestyle='--', x_inline=False, y_inline=False)
gl.top_labels = False
gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.rotate_labels = False
gl.xlocator = mticker.FixedLocator(np.arange(117, 128, 2))
gl.ylocator = mticker.FixedLocator(np.arange(20, 28, 1))
gl.xlabel_style = {'size': 12}
gl.ylabel_style = {'size': 12}
return ax
# Color settings
rgb = ([237, 237, 237], [209, 209, 209], [173, 173, 173], [131, 131, 131],
[93, 93, 93], [151, 198, 223], [111, 176, 214], [49, 129, 189],
[26, 104, 174], [8, 79, 153], [62, 168, 91], [110, 193, 115],
[154, 214, 149], [192, 230, 185], [223, 242, 217], [255, 255, 164],
[255, 243, 0], [255, 183, 0], [255, 123, 0], [255, 62, 0],
[255, 2, 0], [196, 0, 0], [136, 0, 0])
colors = np.array(rgb)/255.
cmap = ListedColormap(colors)
rain_levels = [0.1, 1, 2, 5, 7.5, 10, 13, 16, 20, 25, 30, 35, 40,
50, 60, 70, 80, 90, 100, 125, 150, 175, 200, 250]
rain_levels = [value * 5 for value in rain_levels]
# Read data
ncfile = Dataset('/data/wrfout_d01_2008-07-27_21_00_00')
lon = np.array(ncfile['XLONG'])[0,:,:]
lat = np.array(ncfile['XLAT'])[0,:,:]
R = (to_np(getvar(ncfile, "RAINC")) +
to_np(getvar(ncfile, "RAINNC")) +
to_np(getvar(ncfile, "RAINSH")))
slp = getvar(ncfile, "slp")
smooth_slp = smooth2d(slp, 3, cenweight=4)
# Plot rainfall
fig1, ax1 = plt.subplots(figsize=(10, 8), subplot_kw={'projection': ccrs.PlateCarree()})
cf1 = ax1.contourf(lon, lat, R, levels=rain_levels, cmap=cmap,
transform=ccrs.PlateCarree(), extend='max')
ax1 = setup_map(ax1, [117, 128, 20, 28])
cbar1 = fig1.colorbar(cf1, ax=ax1, ticks=rain_levels[::3], shrink=0.65)
cbar1.set_label('Rainfall (mm)', size=12)
plt.title('Accumulated Rainfall', fontsize=14, pad=20)
plt.savefig('rainfall_plot.png', dpi=300, bbox_inches='tight')
# Plot SLP
fig2, ax2 = plt.subplots(figsize=(10, 8), subplot_kw={'projection': ccrs.PlateCarree()})
cf2 = ax2.contourf(lon, lat, to_np(smooth_slp), levels=10,
transform=ccrs.PlateCarree())
ax2 = setup_map(ax2, [117, 128, 20, 28])
cbar2 = fig2.colorbar(cf2, ax=ax2, shrink=0.65)
cbar2.set_label('Sea Level Pressure (hPa)', size=12)
plt.title('Sea Level Pressure', fontsize=14, pad=20)
plt.savefig('slp_plot.png', dpi=300, bbox_inches='tight')
plt.close('all')降水对比
地形减半降水： 
正常地形降水: 
海平面气压对比
地形减半气压: 
正常地形气压: 
台风路径对比如下： 
五、实验结论
通过比较不同地形条件下的模拟结果，可以分析：
1.台风路径的变化： 降低地形高度后，台风的路径发生了明显改变。与原始地形条件下的路径相比，减半地形条件下的台风路径整体偏南，并且在某些时段出现了打转现象，这表明地形对台风的引导气流和移动方向具有显著影响。
2.台风强度的影响： 海平面气压（SLP）的对比图显示，地形高度的改变对台风的强度也产生了影响。在地形减半的模拟中，台风的最低海平面气压更低，这暗示了地形对台风结构和强度的调制作用。
3.降水分布的差异： 累积降水量的对比图揭示了地形对降水分布和强度的重要影响。在地形减半的模拟中，降水大值区集中在台湾省南部，且降水极值相对于正常地形更大。而在原始地形区域，降水分布较为均匀。
综上所述，地形在台风的路径、强度和降水分布中扮演着关键角色。本实验通过地形敏感性分析，验证了WRF模式在模拟台风与地形相互作用方面的能力，并强调了地形数据精度在台风预报中的重要性。