本文主要介紹Takaya and Nakamura (2001)推導(dǎo)的波活動(dòng)通量（wave activity flux, WAF）的Python實(shí)現(xiàn)。

主要參考了一些大佬完成的WAF水平分量的Python腳本；而由于垂直分量的使用率較低，幾乎找不到現(xiàn)成的Python腳本，故只能根據(jù)日本學(xué)者提供的NCL腳本來改寫。關(guān)于WAF的定義和用途我們不做過多介紹，直接上公式：

計(jì)算WAF的氣象物理量有溫度氣候態(tài)Tc（K）、緯向風(fēng)氣候態(tài)Uc（m/s）、經(jīng)向風(fēng)氣候態(tài)Vc（m/s）、位勢(shì)氣候態(tài)GEOc（m2/s2）和位勢(shì)擾動(dòng)場(chǎng)GEOa。實(shí)際應(yīng)用中位勢(shì)擾動(dòng)場(chǎng)可以用位勢(shì)擾動(dòng)的合成場(chǎng)或回歸場(chǎng)等代替，以計(jì)算合成場(chǎng)或回歸場(chǎng)的WAF。

我們?cè)O(shè)計(jì)函數(shù)TN_flux(Tc,Uc,Vc,GEOc,GEOa)，其功能為：所有輸入變量均為shape相同的(level,lat,lon)三維DataArray。如果level維的長(zhǎng)度是1，不能計(jì)算垂直分量，只輸出水平分量；如果level維的長(zhǎng)度大于1，則輸出三維分量。

由于水平分量相對(duì)簡(jiǎn)單，很多公眾號(hào)的推文給出了解析和腳本，而且文章最后也提供了完整腳本，所以我們不再贅述，這里只詳細(xì)剖析垂直分量。

f0為科氏參數(shù)，公式為f0=2Ωsinφ，其中Ω為地球自轉(zhuǎn)角速度7.292e-5 s^-1，主要代碼為

N^2是Brunt-Vaisala頻率，我們直接看NCL的主要代碼

可以發(fā)現(xiàn)代碼做了一定的簡(jiǎn)化，大氣標(biāo)高設(shè)置為8000m，氣壓標(biāo)高公式為

dthetadz是位溫對(duì)高度的差分，高度可以根據(jù)氣壓推算，也是近似值；gc為氣體常數(shù)，下面直接上python主要代碼

后面擾動(dòng)流函數(shù)對(duì)高度的差分舉一反三，也不再贅述。最后給出完整的腳本，值得注意的是計(jì)算時(shí)保持維度的統(tǒng)一可以避免出錯(cuò)，這里隨時(shí)保持(level,lat,lon)的維度順序，如果是一維、二維數(shù)據(jù)可以用numpy.newaxis擴(kuò)展到三維：

PS：感謝@玄能改命-氪必救非同學(xué)訂正了腳本中的一些錯(cuò)誤。

參考文獻(xiàn)和鏈接

Takaya, K., and H. Nakamura, 2001. A formulation of a phase-independent wave-activity flux for stationary and migratory quasigeostrophic eddies on a zonally varying basic flow. Journal of the Atmospheric Sciences, 58(6): 608-627.

?

日本學(xué)者提供的GrADS、Fortran和 NCL腳本

http://www.atmos.rcast.u-tokyo.ac.jp/nishii/programs/index.html

?

氣象水文科研貓：Python下載ERA5數(shù)據(jù)并計(jì)算T-N Flux波作用通量

https://mp.weixin.qq.com/s/Q1nl-n1tXUq3Tw1PvQPbrQ