大氣動力學(xué)中通常用波作用通量來診斷 Rossby波的傳播睹耐。常用的三種波作用通量分別為局地E-P 通量，Plumb 波作用通量和T-N 波作用通量。局地E-P 通量可以診斷一段時(shí)間內(nèi)天氣尺度瞬變波對定長波的調(diào)控作用徘禁，但無法直接表現(xiàn)Rossby長波的時(shí)間演變過程疾渴，然而T-N 波或Plumb 波作用通量可以 千贯，因此需要診斷一次天氣過程的波活動演變特征時(shí)，通常選擇T-N 波或Plumb 波作用通量搞坝。
Plumb波作用通量與T-N波作用通量的優(yōu)劣
Plumb波作用通量
Plumb(1985)使用小振幅定常波在緯向均勻基本流中傳播時(shí)的守恒關(guān)系搔谴，給出了定常Rossby波的三維波活動通量，代表波能量的傳播方向桩撮。Plumb 波作用通量提出后被廣泛使用敦第，能有效分析定常Rossby波的三維傳播特性。Plumb 波作用通量的緯向分量較大而經(jīng)向分量較小店量，適用于振幅較小的緯向均勻的西風(fēng)帶Rossby長波的診斷芜果。
T-N波作用通量
Takaya 和 Nakamura( 2001 ) 為了更好地診斷真實(shí)大氣中Rossby波的三維傳播特征，將Plumb 波通量進(jìn)一步推廣融师，使其更適用于復(fù)雜的背景氣流右钾，發(fā)展了T-N波作用通量。T-N 波作用通量是對Plumb 波作用通量的改進(jìn)，經(jīng)向分量得以增強(qiáng)舀射，能更好地描述緯向非均勻氣流中的較大振幅的西風(fēng)帶Rossby 長波擾動窘茁，但需要根據(jù)研究目的人為選擇背景場。
數(shù)據(jù)準(zhǔn)備
計(jì)算Plumb通量脆烟，需要準(zhǔn)備所需時(shí)間的位勢高度場山林，UV風(fēng)場；
計(jì)算T-N通量浩淘，需要所需時(shí)間的位勢高度場以及所處季節(jié)(月份)的多年氣候平均場捌朴。
這里我以2014年11月的平均通量場為例(選取這一時(shí)間是由于有正確的圖像作為對比，日本學(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位勢高度場左驾，UV風(fē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ì)算公式：
其中上標(biāo)帶“－”和“ ' ”的變量分別表示緯圈平均和緯向偏差轻猖。φ帆吻，λ，Φ分別表示緯度咙边，經(jīng)度和位勢猜煮，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))
#微分部分办陷，對經(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ì)算垂直分量兩個(gè)原因，一個(gè)是我暫時(shí)沒有正確的例子對比制圈，另一個(gè)是我目前沒有整層大氣的溫度數(shù)據(jù)们童。
T-N波作用通量計(jì)算的Python實(shí)現(xiàn)
Ψ= Φ/ f 是準(zhǔn)地轉(zhuǎn)流函數(shù)相對于氣候場的擾動【校基本流場U = ( U慧库，V ) 表示氣候場，其余參數(shù)與Plumb計(jì)算相同馋嗜。
因此可以發(fā)現(xiàn)齐板，T-N的計(jì)算結(jié)果與背景場的選擇有重要的關(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ì)算括號外的參數(shù)，a^2可以從括號內(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ì)算應(yīng)該沒有什么錯(cuò)誤。但從圖來看莺债，T-N波似乎更加平滑和流暢滋觉，plumb波輻合輻散的特征明顯一點(diǎn)，當(dāng)然了九府，具體的分析還是要結(jié)合天氣事件各個(gè)槽脊系統(tǒng)或者波列的發(fā)展傳播來分析椎瘟。
致謝：感謝dd_atmo11的bug反饋，代碼已更正侄旬。
