python--由wrfouput的数据计算位势涡度，并插值到指定压力层

近日，需要对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])

这样就得到一个文件的位涡啦