近日，需要对wrf模式输出的数据进行计算位涡，并绘图分析。发现模式本身输出的数据中虽然不包含位涡，但wrf-python 提供了函数可以通过其他变量计算得到位涡。顺便记录一下计算的过程以及将位涡插值到压力层的过程
位涡的计算主要通过wrf-python这个库中的wrf.pvo这个函数，官网链接如下：
wrf.pvo(ustag, vstag, theta, pres, msfu, msfv, msfm, cor, dx, dy, meta=True)
计算位涡需要用到的变量包括：
- 纬向风速U
- 经向风速V
- 位势温度T
- 压力P
- 几个其他参数：msfu、msfv、msfm、cor
- x格点之间的距离：dx
- y格点之间的距离：dy
函数计算返回的就是位涡，单位为（puv）1 PVU = 1.0 x 10^(-6) m^2 s^(-1) K kg^(-1)
比较了一下ncl中的计算函数wrf_pvo（）中计算过程是需要将：位涡+300 的，但是wrf-python中的没有体现，按道理应该是以ncl中的为准。
后来发现，通过getvar("T")得到的位温，它是perturbation potential temperature theta-t0，加上300才是 total potential temperature。
下面给出定义的函数，直接输入文件路径，即可得到计算得到的位涡：
def cal_interp(file): ######################################################################## ####通过函数计算位涡 ####################################################################### ncfile = Dataset(filelist[0]) U = getvar(ncfile, "U") V = getvar(ncfile, "V") Theta = getvar(ncfile, "T") P = getvar(ncfile, "P") PB = getvar(ncfile, "PB") MSFU = getvar(ncfile, "MAPFAC_U") MSFV = getvar(ncfile, "MAPFAC_V") MSFM = getvar(ncfile, "MAPFAC_M") COR = getvar(ncfile, "F") DX = ncfile.DX DY = ncfile.DY THETA = Theta + 300 P = P + PB pv = pvo(U, V, THETA, P, MSFU, MSFV, MSFM, COR, DX, DY, 0) ######################################################################## ####下面将得到的位涡插值到等压面上，选取常用的27层压力层进行插值 ####################################################################### pre = getvar(ncfile, "pressure") level =np.array( [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] ) level=level[::-1] plevs =level*units.hPa pv_interp = np.array(interplevel(pv, pre, plevs)) potential_vorticity=pv_interp*units['K m**2 kg**-1 s**-1'] return potential_vorticity
使用函数：
path=r'/wrfout/' filelist = glob.glob(path+'*') filelist.sort() pv = cal_interp(filelist[0])
这样就得到一个文件的位涡啦
