近日，需要對wrf模式輸出的數據進行計算位渦，並繪圖分析。發現模式本身輸出的數據中雖然不包含位渦，但wrf-python 提供了函數可以通過其他變量計算得到位渦。順便記錄一下計算的過程以及將位渦插值到壓力層的過程

位渦的計算主要通過wrf-python這個庫中的wrf.pvo這個函數，官網鏈接如下：

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])

這樣就得到一個文件的位渦啦