大氣動(dòng)力學(xué)中通常用波作用通量來診斷 Rossby波的傳播。常用的三種波作用通量分別為局地E-P 通量，Plumb 波作用通量和T-N 波作用通量。局地E-P 通量可以診斷一段時(shí)間內(nèi)天氣尺度瞬變波對(duì)定長(zhǎng)波的調(diào)控作用，但無法直接表現(xiàn)Rossby長(zhǎng)波的時(shí)間演變過程，然而T-N 波或Plumb 波作用通量可以 ，因此需要診斷一次天氣過程的波活動(dòng)演變特征時(shí)，通常選擇T-N 波或Plumb 波作用通量。

2014年11月300hPa兩種波作用通量

Plumb波作用通量與T-N波作用通量的優(yōu)劣

Plumb波作用通量

Plumb(1985)使用小振幅定常波在緯向均勻基本流中傳播時(shí)的守恒關(guān)系，給出了定常Rossby波的三維波活動(dòng)通量，代表波能量的傳播方向。Plumb 波作用通量提出后被廣泛使用，能有效分析定常Rossby波的三維傳播特性。Plumb 波作用通量的緯向分量較大而經(jīng)向分量較小，適用于振幅較小的緯向均勻的西風(fēng)帶Rossby長(zhǎng)波的診斷。

T-N波作用通量

Takaya 和 Nakamura( 2001 ) 為了更好地診斷真實(shí)大氣中Rossby波的三維傳播特征，將Plumb 波通量進(jìn)一步推廣，使其更適用于復(fù)雜的背景氣流，發(fā)展了T-N波作用通量。T-N 波作用通量是對(duì)Plumb 波作用通量的改進(jìn)，經(jīng)向分量得以增強(qiáng)，能更好地描述緯向非均勻氣流中的較大振幅的西風(fēng)帶Rossby 長(zhǎng)波擾動(dòng)，但需要根據(jù)研究目的人為選擇背景場(chǎng)。

數(shù)據(jù)準(zhǔn)備

計(jì)算Plumb通量，需要準(zhǔn)備所需時(shí)間的位勢(shì)高度場(chǎng)，UV風(fēng)場(chǎng)；
計(jì)算T-N通量，需要所需時(shí)間的位勢(shì)高度場(chǎng)以及所處季節(jié)(月份)的多年氣候平均場(chǎng)。
這里我以2014年11月的平均通量場(chǎng)為例(選取這一時(shí)間是由于有正確的圖像作為對(duì)比，日本學(xué)者提供了相應(yīng)的Grads及fortran計(jì)算腳本
http://www.atmos.rcast.u-tokyo.ac.jp/nishii/programs/index.html)。

#以下為數(shù)據(jù)讀取部分，最終所得z300，u300，v300為2014年11月的300hPa位勢(shì)高度場(chǎng)，UV風(fēng)場(chǎng)的月平均；
#z_tmp，u_tmp，v_tmp為1979-2018年11月的氣候態(tài)；
#所有數(shù)據(jù)來自ECMWF的ERA-Interim數(shù)據(jù)集
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波作用通量計(jì)算的Python實(shí)現(xiàn)

首先先給出Plumb波作用通量的計(jì)算公式：

Plumb波作用通量的計(jì)算公式

其中上標(biāo)帶“－”和“ ' ”的變量分別表示緯圈平均和緯向偏差。φ，λ，Φ分別表示緯度，經(jīng)度和位勢(shì)，f = 2Ωsinφ表示科氏參數(shù)，a，Ω 分別表示地球半徑和地球自轉(zhuǎn)速率。p0 = 1000 hPa 。目前暫時(shí)不涉及垂直分量。
計(jì)算難點(diǎn)在于偏微分部分，思路就在于將微分變?yōu)椴罘钟?jì)算，Numpy提供了關(guān)于差分的方法有兩個(gè)，一個(gè)是使用np.diff()函數(shù)計(jì)算前后差，另一個(gè)方法是使用np.gradient()計(jì)算相鄰梯度，然后除以格距。這里我使用的是np.gradient()。

a=6400000 #地球半徑
omega=7.292e-5 #自轉(zhuǎn)角速度
lev = 300/1000#p/p0
#要把經(jīng)緯度轉(zhuǎn)換成角度量，所以做（*np.pi/180.0）處理
#因?yàn)樽罱K要計(jì)算Fx,Fy，所以統(tǒng)一數(shù)組shape，使用.reshape((1,-1))或(-1,1)處理
#只有經(jīng)度維的使用((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))
#計(jì)算緯偏值
ua = u300 - u300.mean((1)).reshape((-1,1))
va = v300 - v300.mean((1)).reshape((-1,1))
za = z300 - z300.mean(1).reshape((-1,1))
#微分部分，對(duì)經(jīng)度求微分，因此是axis=1，dlon是之前求的經(jīng)度格距
duzdlon = np.gradient(ua * za ,axis = 1)/dlon
dvzdlon = np.gradient(va * za ,axis = 1)/dlon
#根據(jù)公式，把各個(gè)部件組合起來，計(jì)算完畢
fx = p_p0*coslat*(va*va - dvzdlon/(sin2lat*2*omega*a))
fy = p_p0*coslat*((-1)*ua*va + duzdlon/(sin2lat*2*omega*a))

我沒計(jì)算垂直分量?jī)蓚€(gè)原因，一個(gè)是我暫時(shí)沒有正確的例子對(duì)比，另一個(gè)是我目前沒有整層大氣的溫度數(shù)據(jù)。

T-N波作用通量計(jì)算的Python實(shí)現(xiàn)

T-N

Ψ= Φ/ f 是準(zhǔn)地轉(zhuǎn)流函數(shù)相對(duì)于氣候場(chǎng)的擾動(dòng)?；玖鲌?chǎng)U = ( U，V ) 表示氣候場(chǎng)，其余參數(shù)與Plumb計(jì)算相同。
因此可以發(fā)現(xiàn)，T-N的計(jì)算結(jié)果與背景場(chǎng)的選擇有重要的關(guān)系。
T-N的計(jì)算難度在于二階微分的計(jì)算，其實(shí)就是差分的過程計(jì)算兩次。

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

#計(jì)算科氏力
f=np.array(2*omega*np.sin(lat*np.pi/180.0)).reshape((-1,1)) 
#計(jì)算|U|
wind = np.sqrt(u**2+v**2)
#計(jì)算括號(hào)外的參數(shù)，a^2可以從括號(hào)內(nèi)提出
c=(lev)*coslat/(2*a*a*wind)
#Φ`
za = z300 - z_tmp.mean(1).reshape((-1,1))
#Ψ`
streamf = g*za/f
#計(jì)算各個(gè)部件，難度在于二階導(dǎo)，變量的名字應(yīng)該可以很容易看出我是在計(jì)算哪部分
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
#計(jì)算兩個(gè)分量
fx = c * ((u/coslat/coslat)*x_tmp+v*y_tmp/coslat)
fy = c * ((u/coslat)*y_tmp+v*x_tmp)

我計(jì)算得到的

這是日本學(xué)者提供的

兩者對(duì)比來看我的計(jì)算應(yīng)該沒有什么錯(cuò)誤。但從圖來看，T-N波似乎更加平滑和流暢，plumb波輻合輻散的特征明顯一點(diǎn)，當(dāng)然了，具體的分析還是要結(jié)合天氣事件各個(gè)槽脊系統(tǒng)或者波列的發(fā)展傳播來分析。

致謝：感謝dd_atmo11的bug反饋，代碼已更正。