如何用wrfout计算水汽通量散度 - 成就云开发者社区

前言

本文旨在实现WRFOUT的单层水汽通量散度和整层水汽通量散度计算方法。WRF（Weather Research and Forecasting）模式是一种广泛应用于天气和气候预测研究的数值模式。水汽通量散度在天气和气候研究中具有重要作用。本项目将针对WRF模式的输出数据（WRFOUT）进行处理和分析，实现单层水汽通量散度和整层水汽通量散度的计算。

在实现该功能的过程中，下面将详细介绍所采用的公式原理，并给出相应的代码示例和使用说明。同时会对计算结果进行可视化展示，以便更好地理解和分析水汽通量散度的空间分布和变化规律。

概念简介

水汽通量散度是衡量水汽输送量变化的一个指标。

水汽通量散度表示单位时间内和单位面积上的水汽通量变化率。它可以反映水汽是否聚集或分散,更能准确地判断是否有利于对流天气的形成。

水汽通量散度的负值表示水汽聚集,这有利于对流天气的形成;反之,正值表示水汽分散,不利于对流天气。

水汽通量散度公式

本文计算部分参考了

https://blog.csdn.net/weixin_44237337/article/details/122601116

单层与整层的概念可以阅读

https://blog.csdn.net/wdbhysszjswn/article/details/129044910

单层水汽实现

1.导入库

代码语言：javascript

复制

代码语言：javascript

复制

from wrf import uvmet, to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords
import numpy as np
from netCDF4 import Dataset
import xarray as xr
from metpy.units import units
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from matplotlib.colors import from_levels_and_colors
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import cmaps
from glob import glob
import metpy.calc as mpcalc
import metpy.constants as constants

代码语言：javascript

复制

代码语言：javascript

复制

Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1

2. 读取数据

代码语言：javascript

复制

代码语言：javascript

复制

# 用 netCDF4 包读取 WRF 模拟数据
wrf_file = Dataset('/home/mw/input/wrfout3385/wrfout_d02_2022-07-14_0600.nc')
# 提取要素
lon = getvar(wrf_file, 'lon', timeidx=0)
lat = getvar(wrf_file, 'lat', timeidx=0)
u = getvar(wrf_file, 'ua', timeidx=0)
v = getvar(wrf_file, 'va', timeidx=0)
q = getvar(wrf_file, 'QVAPOR', timeidx=0)*1000
# 提取WRF模拟的经纬度数组
lats, lons = latlon_coords(u)
提取WRF模拟的地图投影
wrf_proj = get_cartopy(u)

代码语言：javascript

复制

3. 散度计算过程

代码语言：javascript

复制

代码语言：javascript

复制

p = getvar(wrf_file, 'pressure', timeidx=0)

u850 = interplevel(u, p, 850)

v850 = interplevel(v, p, 850)

q850 = interplevel(q, p, 850) ** units('g/kg')

z = getvar(wrf_file, 'z', units="dm", timeidx=0)

z850 = interplevel(z, p, 850)

dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
计算水汽通量散度
q_flux_divergence = mpcalc.divergence(

to_np(

u850 *

q850/constants.g),

to_np(

v850 *

q850/constants.g),

dx=to_np(dx),

dy=to_np(dy))

代码语言：javascript

复制

4. 绘图

整层水汽

代码语言：javascript

复制

代码语言：javascript

复制

代码语言：javascript

复制

u_interp = interplevel(u, p, levels)

v_interp = interplevel(v, p, levels)

q_interp = interplevel(q, p, levels) ** units('g/kg')

lev = u_interp.level
计算水汽通量散度
qu = u_interp *q_interp/constants.g

qv = v_interp *q_interp/constants.g

q_flux_divergence_all = np.zeros((lev.shape[0],lat.shape[0],lon.shape[0]))

for i in range(lev.shape[0]):

q_flux_divergence_all[i] = mpcalc.divergence(u = to_np(qu[i]),v = to_np(qv[i]),dx = to_np(dx) ,dy = to_np(dy),x_dim=-1, y_dim=-2)
将 q_flux_divergence_all 中的 NaN 值替换为 0
q_flux_divergence_all = np.nan_to_num(q_flux_divergence_all, nan=0)

total_div_qv = np.trapz(q_flux_divergence_all,lev,axis=0)

total_div_qv[2]

最后

metpy要注意的点挺多的，什么单位，维度，格式都要对齐才能好好运行，大家只能多测试。

完整代码与文件可后台回复“水汽通量散度”获取链接。

封面由ai绘制

/imagine prompt: color photo of water vapor flux, depicting the movement and exchange of moisture in the atmosphere —c 10 —ar 2:3