气象数据处理：T-N三维波活动通量的Python实现

本文主要介绍Takaya and Nakamura (2001)推导的波活动通量（wave activity flux, WAF）的Python实现。

主要参考了一些大佬完成的WAF水平分量的Python脚本；而由于垂直分量的使用率较低，几乎找不到现成的Python脚本，故只能根据日本学者提供的NCL脚本来改写。关于WAF的定义和用途我们不做过多介绍，直接上公式：

计算WAF的气象物理量有温度气候态Tc（K）、纬向风气候态Uc（m/s）、经向风气候态Vc（m/s）、位势气候态GEOc（m²/s²）和位势扰动场GEOa。实际应用中位势扰动场可以用位势扰动的合成场或回归场等代替，以计算合成场或回归场的WAF。

我们设计函数TN_flux(Tc,Uc,Vc,GEOc,GEOa)，其功能为：所有输入变量均为shape相同的(level,lat,lon)三维DataArray。如果level维的长度是1，不能计算垂直分量，只输出水平分量；如果level维的长度大于1，则输出三维分量。

由于水平分量相对简单，很多公众号的推文给出了解析和脚本，而且文章最后也提供了完整脚本，所以我们不再赘述，这里只详细剖析垂直分量。

f0为科氏参数，公式为f0=2Ωsinφ，其中Ω为地球自转角速度7.292e-5 s^-1，主要代码为

N^2是Brunt-Vaisala频率，我们直接看NCL的主要代码

可以发现代码做了一定的简化，大气标高设置为8000m，气压标高公式为

dthetadz是位温对高度的差分，高度可以根据气压推算，也是近似值；gc为气体常数，下面直接上python主要代码

后面扰动流函数对高度的差分举一反三，也不再赘述。最后给出完整的脚本，值得注意的是计算时保持维度的统一可以避免出错，这里随时保持(level,lat,lon)的维度顺序，如果是一维、二维数据可以用numpy.newaxis扩展到三维：

PS：感谢@玄能改命-氪必救非同学订正了脚本中的一些错误。

参考文献和链接

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.

 

日本学者提供的GrADS、Fortran和 NCL脚本

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

 

气象水文科研猫：Python下载ERA5数据并计算T-N Flux波作用通量

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