大气动力学中通常用波作用通量来诊断 Rossby波的传播。常用的三种波作用通量分别为局地E-P 通量，Plumb 波作用通量和T-N 波作用通量。局地E-P 通量可以诊断一段时间内天气尺度瞬变波对定长波的调控作用，但无法直接表现Rossby长波的时间演变过程，然而T-N 波或Plumb 波作用通量可以 ，因此需要诊断一次天气过程的波活动演变特征时，通常选择T-N 波或Plumb 波作用通量。

2014年11月300hPa两种波作用通量

Plumb波作用通量与T-N波作用通量的优劣

Plumb波作用通量

Plumb(1985)使用小振幅定常波在纬向均匀基本流中传播时的守恒关系，给出了定常Rossby波的三维波活动通量，代表波能量的传播方向。Plumb 波作用通量提出后被广泛使用，能有效分析定常Rossby波的三维传播特性。Plumb 波作用通量的纬向分量较大而经向分量较小，适用于振幅较小的纬向均匀的西风带Rossby长波的诊断。

T-N波作用通量

Takaya 和 Nakamura( 2001 ) 为了更好地诊断真实大气中Rossby波的三维传播特征，将Plumb 波通量进一步推广，使其更适用于复杂的背景气流，发展了T-N波作用通量。T-N 波作用通量是对Plumb 波作用通量的改进，经向分量得以增强，能更好地描述纬向非均匀气流中的较大振幅的西风带Rossby 长波扰动，但需要根据研究目的人为选择背景场。

数据准备

计算Plumb通量，需要准备所需时间的位势高度场，UV风场；
计算T-N通量，需要所需时间的位势高度场以及所处季节(月份)的多年气候平均场。
这里我以2014年11月的平均通量场为例(选取这一时间是由于有正确的图像作为对比，日本学者提供了相应的Grads及fortran计算脚本
http://www.atmos.rcast.u-tokyo.ac.jp/nishii/programs/index.html)。

#以下为数据读取部分，最终所得z300，u300，v300为2014年11月的300hPa位势高度场，UV风场的月平均；
#z_tmp，u_tmp，v_tmp为1979-2018年11月的气候态；
#所有数据来自ECMWF的ERA-Interim数据集
f_z300 = xr.open_dataset("z.nc")
z300 = np.array(f_z300['z'].loc[:,300,:,:].loc['2014-11-01'])
z300 = z300/9.8
z_tmp = np.array(f_z300['z'].loc[:,300,:,:].loc[f_z300.time.dt.month.isin([11])])
z_tmp = z_tmp.mean((0))/9.8

f_u300 = xr.open_dataset("u.nc")
u300 = np.array(f_u300['u'].loc[:,300,:,:].loc['2014-11-01'])
u_tmp = np.array(f_u300['u'].loc[:,300,:,:].loc[f_u300.time.dt.month.isin([11])])

f_v300 = xr.open_dataset("v.nc")
v300 = np.array(f_v300['v'].loc[:,300,:,:].loc['2014-11-01'])
v_tmp = np.array(f_v300['v'].loc[:,300,:,:].loc[f_v300.time.dt.month.isin([11])])

lat = f_z300['latitude']
lon = f_z300['longitude']

Plumb波作用通量计算的Python实现

首先先给出Plumb波作用通量的计算公式：

Plumb波作用通量的计算公式

其中上标带“－”和“ ' ”的变量分别表示纬圈平均和纬向偏差。φ，λ，Φ分别表示纬度，经度和位势，f = 2Ωsinφ表示科氏参数，a，Ω 分别表示地球半径和地球自转速率。p0 = 1000 hPa 。目前暂时不涉及垂直分量。
计算难点在于偏微分部分，思路就在于将微分变为差分计算，Numpy提供了关于差分的方法有两个，一个是使用np.diff()函数计算前后差，另一个方法是使用np.gradient()计算相邻梯度，然后除以格距。这里我使用的是np.gradient()。

a=6400000 #地球半径
omega=7.292e-5 #自转角速度
lev = 300/1000#p/p0
#要把经纬度转换成角度量，所以做（*np.pi/180.0）处理
#因为最终要计算Fx,Fy，所以统一数组shape，使用.reshape((1,-1))或(-1,1)处理
#只有经度维的使用((1,-1))，只有纬度维的使用((-1,1))
dlon=(np.gradient(lon)*np.pi/180.0).reshape((1,-1))
coslat = (np.cos(np.array(lat)*np.pi/180)).reshape((-1,1))
sin2lat = (np.sin(np.array(lat)*2*np.pi/180)).reshape((-1,1))
#计算纬偏值
ua = u300 - u300.mean((1)).reshape((-1,1))
va = v300 - v300.mean((1)).reshape((-1,1))
za = z300 - z300.mean(1).reshape((-1,1))
#微分部分，对经度求微分，因此是axis=1，dlon是之前求的经度格距
duzdlon = np.gradient(ua * za ,axis = 1)/dlon
dvzdlon = np.gradient(va * za ,axis = 1)/dlon
#根据公式，把各个部件组合起来，计算完毕
fx = p_p0*coslat*(va*va - dvzdlon/(sin2lat*2*omega*a))
fy = p_p0*coslat*((-1)*ua*va + duzdlon/(sin2lat*2*omega*a))

我没计算垂直分量两个原因，一个是我暂时没有正确的例子对比，另一个是我目前没有整层大气的温度数据。

T-N波作用通量计算的Python实现

T-N

Ψ= Φ/ f 是准地转流函数相对于气候场的扰动。基本流场U = ( U，V ) 表示气候场，其余参数与Plumb计算相同。
因此可以发现，T-N的计算结果与背景场的选择有重要的关系。
T-N的计算难度在于二阶微分的计算，其实就是差分的过程计算两次。

dlon=(np.gradient(lon)*np.pi/180.0).reshape((1,-1))
dlat=(np.gradient(lat)*np.pi/180.0).reshape((-1,1))
coslat = (np.cos(np.array(lat)*np.pi/180)).reshape((-1,1))
sinlat = (np.sin(np.array(lat)*np.pi/180)).reshape((-1,1))

#计算科氏力
f=np.array(2*omega*np.sin(lat*np.pi/180.0)).reshape((-1,1)) 
#计算|U|
wind = np.sqrt(u**2+v**2)
#计算括号外的参数，a^2可以从括号内提出
c=(lev)*coslat/(2*a*a*wind)
#Φ`
za = z300 - z_tmp.mean(1).reshape((-1,1))
#Ψ`
streamf = g*za/f
#计算各个部件，难度在于二阶导，变量的名字应该可以很容易看出我是在计算哪部分
dzdlon = np.gradient(streamf,axis = 1)/dlon
ddzdlonlon = np.gradient(dzdlon,axis = 1)/dlon
dzdlat = np.gradient(streamf,axis = 0)/dlat
ddzdlatlat = np.gradient(dzdlat,axis = 0)/dlat
ddzdlatlon = np.gradient(dzdlat,axis = 1)/dlon
#这是X,Y分量共有的部分
x_tmp = dzdlon*dzdlon-streamf*ddzdlonlon
y_tmp = dzdlon*dzdlat-streamf*ddzdlatlon
#计算两个分量
fx = c * ((u/coslat/coslat)*x_tmp+v*y_tmp/coslat)
fy = c * ((u/coslat)*y_tmp+v*x_tmp)

我计算得到的 这是日本学者提供的

两者对比来看我的计算应该没有什么错误。但从图来看，T-N波似乎更加平滑和流畅，plumb波辐合辐散的特征明显一点，当然了，具体的分析还是要结合天气事件各个槽脊系统或者波列的发展传播来分析。

致谢：感谢dd_atmo11的bug反馈，代码已更正。