近日,需要對wrfModel output data to calculate potential vorticity,並繪圖分析.It is found that the data output by the mode itself does not contain potential vortices,但wrf-python Provides a function to calculate potential vorticity with other variables.By the way, record the process of calculation and the process of interpolating the potential vorticity into the pressure layer

The calculation of potential vorticity is mainly passedwrf-python這個庫中的wrf.pvo這個函數,官網鏈接如下：

wrf-pvo

wrf.pvo(ustag, vstag, theta, pres, msfu, msfv, msfm, cor, dx, dy, meta=True)

The variables needed to calculate the potential vorticity include：

• 緯向風速U
• 經向風速V
• Potential temperatureT
• 壓力P
• several other parameters：msfu、msfv、msfm、cor
• xdistance between grid points：dx
• ydistance between grid points：dy
  What the function calculation returns is the potential vorticity,單位為（puv）1 PVU = 1.0 x 10^(-6) m^2 s^(-1) K kg^(-1)

比較了一下ncl中的計算函數wrf_pvo（）The calculation process is required to be：potential vorticity+300 的,但是wrf-pythonis not reflected in,Logically it should bencl中的為准.

後來發現,通過getvar("T")obtained bit temperature,它是perturbation potential temperature theta-t0,加上300才是 total potential temperature.

The defined functions are given below,Enter the file path directly,The calculated potential vorticity can be obtained：

def cal_interp(file):
########################################################################
####Calculate potential vorticity by function
#######################################################################
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)
########################################################################
####The resulting potential vorticity is then interpolated onto the isobaric surface,Select Common27The layer pressure layer is interpolated
#######################################################################
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])

In this way, the potential vortex of a file is obtained