python--测试使用不同的方式计算位涡平流项的差异

python计算位涡平流项

主要为了测试计算位涡平流项的准确性，这里使用了两种方法来计算位涡的平流项。

• 方法1：使用metpy函数
• 方法2：编写代码自己计算

一般使用metpy函数的计算结果，已经是等号右边的那一项了。也就是说相当于下图所示：

展开来就是：

$-$($u$$\frac{\partial pv}{\partial x}$+$v$$\frac{\partial pv}{\partial y}$)

使用metpy的函数时，需要注意的点：变量带上单位。

具体的代码如下所示：

import numpy as np
import pandas as pd
import xarray as xr
import proplot as pplt
import glob
from matplotlib.colors import ListedColormap 
from wrf import getvar,pvo,interplevel,destagger,latlon_coords
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from netCDF4 import Dataset
import matplotlib as mpl   
from metpy.units import units
import metpy.constants as constants
import metpy.calc as mpcalc
import cmaps

f_w   = r'D:/wind.nc'
f_pv  = r'D:/pv.nc'

def  cal_dxdy(file):
    ncfile = Dataset(file)
    P = getvar(ncfile, "pressure")
    lats, lons = latlon_coords(P)
    lon = lons[0]
    lon[lon<=0]=lon[lon<=0]+360
    lat = lats[:,0]
    dx, dy = mpcalc.lat_lon_grid_deltas(lon.data, lat.data)
    return lon,lat,dx,dy
path  = r'J:/wrfout/'
filelist = glob.glob(path+'wrf*')
filelist.sort()
lon,lat,dx,dy = cal_dxdy(filelist[-1]) 

dw = xr.open_dataset(f_w)
dp = xr.open_dataset(f_pv)

u = dw.u
v = dw.v
pv = dp.pv

u = u.data*units('m/s')
v = v.data*units('m/s')
units.define('PVU = 1e-6 m^2 s^-1 K kg^-1 = pvu')
pv = pv.data*units('pvu').to('m^2 s^-1 K kg^-1')

lon = dw.lon
lat = dw.lat
time= dw.time
leve= dw.level

xlon,ylat=np.meshgrid(lon,lat)
dlony,dlonx=np.gradient(xlon)
dlaty,dlatx=np.gradient(ylat)
pi=3.14159265
re=6.37e6
dx=re*np.cos(ylat*pi/180)*dlonx*pi/180
dy=re*dlaty*pi/180


adv = np.full((time.shape[0],leve.shape[0],lat.shape[0],lon.shape[0]),np.nan)*units('m/s')*units("m^2 s^-1 K kg^-1")/units("m")

for i in range(time.shape[0]):
    
    for j in range(leve.shape[0]):
        
        print(i,j) 
        adv[i,j,:,:] = mpcalc.advection(pv[i,j,:,:],u=u[i,j,:,:],v=v[i,j,:,:],dx=dx,dy=dy,x_dim=-1,y_dim=-2)
       
########################## way 2 ########################################## 
u0 = dw.u.data
v0 = dw.v.data
dpdx = np.gradient(pv,axis=-1)
dpdy = np.gradient(pv,axis=-2)

pvadv = np.full((time.shape[0],leve.shape[0],lat.shape[0],lon.shape[0]),np.nan)

for i in range(time.shape[0]):
    
    for j in range(leve.shape[0]):
        
        print(i,j) 
        pvadv[i,j] = -(u0[i,j]*dpdx[i,j]/dx+v0[i,j]*dpdy[i,j]/dy)
        



f, axs = pplt.subplots( ncols=1, nrows=2,
                       figsize=(8,6),
                      tight=True,
                       proj='cyl',
                       proj_kw={'lon_0': 180}, lonlim=(100, 210), 
                             latlim=(-25, 25),
                       share=4,
                       )


axs[0].contourf(lon[::4],lat[::4],adv[16,6][::4,::4],colorbar='b')
axs[1].contourf(lon[::4],lat[::4],pvadv[16,6][::4,::4],colorbar='b')
axs.format(title=time[16],titleloc='l',
          coast=True, 
            labels=True,
           fontsize=10,
          )

结果对比，更方面的是直接绘制折线图。选取两个同样的点即可