Python编程 | T-N波作用通量水平分量
受“甲方”委托，我写了一个计算T-N波作用通量水平分量的Python脚本。虽然之前我从来没有听说过“T-N波作用通量”这个东西，但是好在公式里每个物理量都还算眼熟，仔细捋顺了计算细节，最终成果还是受到了“甲方”的肯定。
单位
通常来说，有公式有数据，总能通过编程算出一个结果，而结果的正确性或者说程序的正确性怎么保证呢？
单位是其中一个很好判断工具。
T-N波作用通量的公式如下
其中
与水平分量有关的物理量如下表所示
符号 | 意义 | 单位 |
|---|---|---|
φ | 纬度（弧度制） | 1 |
λ | 经度（弧度制） | 1 |
a | 地球半径 | m |
p | 气压与1000hPa之比 | 1 |
U | 气候场基本流场U分量 | m/s |
V | 气候场基本流场V分量 | m/s |
|U| | 气候场基本流场 | m/s |
Φ' | 位势的纬向偏差 | m²/s² |
f | 科氏参数 | 1/s |
ψ' | 准地转流函数相对于气候场的扰动 | m²/s² |
由于我也并没有学过相关的内容，以上是综合几篇文献和教程总结或者反推的。比如，一些文献中没有说是气压与1000hPa之比，但是通过T-N波作用通量单位最终为m²/s²以及其他文献判断应当是如此。
位势与位势高度
位势高度的单位是gpm位势米或者m米，gpm在metpy的单位系统中与m是等价的。
位势的单位是m²/s²，数值上约是位势高度的9.8倍。
其实就是位势高度乘以重力加速度等于位势。容易记反的话看一下单位就知道了。
千万要注意你使用的数据是位势还是位势高度！
用位势高度求位势可以用metpy中的height_to_geopotential函数来实现。
偏导
numpy、xarray、metpy都可以求偏导，我其实更喜欢metpy。之前计算水汽通量、Zwack-Okossi诊断方程时都是使用metpy进行梯度（偏导）、二阶偏导、涡度和拉普拉斯等计算，非常方便，但是T-N波作用通量却并不适合用metpy，因为metpy会“自作主张”把对弧度制经纬度的偏导转换为对水平距离的偏导，所以要用xarray或numpy这种相对来说比较“手动”的方案。
numpy的gradient函数和xarray的differentiate方法本质上是一样的，从文档上来看，differentiate应该是底层调用了gradient。用过几次后觉得differentiate更方便，直接指定dim名称即可计算对应的second order accurate central differences（二阶精度中心差商）即偏导。
脚本
完整的脚本如下，过一段时间想办法封装成函数放到GitHub上。
读取数据并标记单位
import numpy as np
import xarray as xr
import metpy.calc as mpcalc
from metpy.units import units
from metpy.constants import earth_avg_radius
p = 250 * units('hPa') # 也有用300hPa的
mon = 1 # 目标月
time_target = f'2015-{mon:02d}-08' # 目标日期
p0 = 1000 * units('hPa')
hgt_day_path = './2015年/hgt.2015.nc'
hgt_mon_path = './2015年/hgt.mon.mean.nc'
uwnd_day_path = './2015年/uwnd.2015.nc'
uwnd_mon_path = './2015年/uwnd.mon.mean.nc'
vwnd_day_path = './2015年/vwnd.2015.nc'
vwnd_mon_path = './2015年/vwnd.mon.mean.nc'
hgt_day = xr.open_dataset(hgt_day_path)['hgt'] * units('m')
hgt_mon = xr.open_dataset(hgt_mon_path)['hgt'] * units('m')
uwnd_day = xr.open_dataset(uwnd_day_path)['uwnd'] * units('m/s')
uwnd_mon = xr.open_dataset(uwnd_mon_path)['uwnd'] * units('m/s')
vwnd_day = xr.open_dataset(vwnd_day_path)['vwnd'] * units('m/s')
vwnd_mon = xr.open_dataset(vwnd_mon_path)['vwnd'] * units('m/s')
metpy的计算函数是接受xarray的dataarray类型的，并且似乎可以自动识别其中的units属性，所以其实并不一定要手动地乘以单位，我这样写主观上觉得更“稳”，另一方面从代码上更便于理解。
# 位势高度转位势
Φ_day = mpcalc.height_to_geopotential(hgt_day)
Φ_mon = mpcalc.height_to_geopotential(hgt_mon)
# 目标日期和目标气压的位势
Φ = Φ_day.sel(time=time_target, level=p)
# 目标月和目标气压位势的气候态
Φ_climatic = Φ_mon.sel(time=Φ_mon['time.month']==mon, level=p).mean(dim='time')
# 目标月和目标气压基本流场U分量的气候态
u_climatic = uwnd_mon.sel(time=uwnd_mon['time.month']==mon, level=p).mean(dim='time')
# 目标月和目标气压基本流场V分量的气候态
v_climatic = vwnd_mon.sel(time=vwnd_mon['time.month']==mon, level=p).mean(dim='time')
# 经纬度转为弧度制
lon_deg = Φ['lon'].copy()
lat_deg = Φ['lat'].copy()
lon_rad = np.deg2rad(lon_deg) * units('1')
lat_rad = np.deg2rad(lat_deg) * units('1')
# 科氏参数
f = mpcalc.coriolis_parameter(Φ['lat'])
# 目标月和目标气压基本流场的气候态
wind_climatic = mpcalc.wind_speed(u_climatic, v_climatic)
cosφ = np.cos(lat_rad)
# 位势的纬向偏差
Φ_prime = Φ - Φ_climatic.mean(dim='lon')
# 将需要对弧度制经纬度求偏导的量的坐标都换成弧度制经纬度
Φ_prime = Φ_prime.assign_coords({'lon': lon_rad, 'lat': lat_rad})
f = f.assign_coords({'lat': lat_rad})
cosφ = cosφ.assign_coords({'lat': lat_rad})
u_climatic = u_climatic.assign_coords({'lon': lon_rad, 'lat': lat_rad})
v_climatic = v_climatic.assign_coords({'lon': lon_rad, 'lat': lat_rad})
wind_climatic = wind_climatic.assign_coords({'lon': lon_rad, 'lat': lat_rad})
# 准地转流函数相对于气候场的扰动
Ψ_prime = Φ_prime / f
# 一顿偏导猛如虎
dΨ_prime_dλ = Ψ_prime.differentiate('lon')
dΨ_prime_dφ = Ψ_prime.differentiate('lat')
ddΨ_prime_ddλ = dΨ_prime_dλ.differentiate('lon')
ddΨ_prime_ddφ = dΨ_prime_dφ.differentiate('lat')
ddΨ_prime_dλφ = dΨ_prime_dλ.differentiate('lat')
# T-N波作用通量的水平分量公共部分
temp1 = p / p0 * cosφ / (2 * wind_climatic * earth_avg_radius**2)
temp2 = dΨ_prime_dλ * dΨ_prime_dφ - Ψ_prime * ddΨ_prime_dλφ
# T-N波作用通量的水平分量
fx = temp1 * (u_climatic / cosφ**2 * (dΨ_prime_dλ**2 - Ψ_prime * ddΨ_prime_ddλ) + v_climatic / cosφ * temp2)
fy = temp1 * (u_climatic / cosφ * temp2 + v_climatic * (dΨ_prime_dφ**2 - Ψ_prime * ddΨ_prime_ddφ))
# 把弧度制经纬度再换成角度制便于画图
lon = lon_deg.values
lon[lon>180] -= 360
fx = fx.assign_coords({'lon': lon, 'lat': lat_deg}).sortby(['lon', 'lat'])
fy = fy.assign_coords({'lon': lon, 'lat': lat_deg}).sortby(['lon', 'lat'])
查看一下fx可以看到，计算出来的单位确实是m²/s²
简单画个图
