Python3.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
Φ'
位势的纬向偏差m2/s2
f科氏参数1/s
ψ'准地转流函数相对于
气候场的扰动m2/s2由于我也并没有学过相关的内容，以上是综合几篇文献和教程总结或者反推的。比如，一些文献中没有说是气压与1000hPa之比，但是通过T-N波作用通量单位最终为m2/s2以及其他文献判断应当是如此。

位势与位势高度

位势高度的单位是gpm位势米或者m米，gpm在metpy的单位系统中与m是等价的。

位势的单位是m2/s2，数值上约是位势高度的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上。

读取数据并标记单位

1. import numpy as np
   
2. import xarray as xr
   
3. import metpy.calc as mpcalc
   
4. from metpy.units import units
   
5. from metpy.constants import earth_avg_radius
   
6. 
   
7. 
   
8. p = 250 * units('hPa')  # 也有用300hPa的
   
9. mon = 1  # 目标月
   
10. time_target = f'2015-{mon:02d}-08'  # 目标日期
    
11. p0 = 1000 * units('hPa')
    
12. 
    
13. hgt_day_path = './2015年/hgt.2015.nc'
    
14. hgt_mon_path = './2015年/hgt.mon.mean.nc'
    
15. uwnd_day_path = './2015年/uwnd.2015.nc'
    
16. uwnd_mon_path = './2015年/uwnd.mon.mean.nc'
    
17. vwnd_day_path = './2015年/vwnd.2015.nc'
    
18. vwnd_mon_path = './2015年/vwnd.mon.mean.nc'
    
19. hgt_day = xr.open_dataset(hgt_day_path)['hgt'] * units('m')
    
20. hgt_mon = xr.open_dataset(hgt_mon_path)['hgt'] * units('m')
    
21. uwnd_day = xr.open_dataset(uwnd_day_path)['uwnd'] * units('m/s')
    
22. uwnd_mon = xr.open_dataset(uwnd_mon_path)['uwnd'] * units('m/s')
    
23. vwnd_day = xr.open_dataset(vwnd_day_path)['vwnd'] * units('m/s')
    
24. vwnd_mon = xr.open_dataset(vwnd_mon_path)['vwnd'] * units('m/s')

复制代码

metpy的计算函数是接受xarray的dataarray类型的，并且似乎可以自动识别其中的units属性，所以其实并不一定要手动地乘以单位，我这样写主观上觉得更“稳”，另一方面从代码上更便于理解。

1. # 位势高度转位势
   
2. Φ_day = mpcalc.height_to_geopotential(hgt_day)
   
3. Φ_mon = mpcalc.height_to_geopotential(hgt_mon)
   
4. # 目标日期和目标气压的位势
   
5. Φ = Φ_day.sel(time=time_target, level=p)
   
6. # 目标月和目标气压位势的气候态
   
7. Φ_climatic = Φ_mon.sel(time=Φ_mon['time.month']==mon, level=p).mean(dim='time')
   
8. # 目标月和目标气压基本流场U分量的气候态
   
9. u_climatic = uwnd_mon.sel(time=uwnd_mon['time.month']==mon, level=p).mean(dim='time')
   
10. # 目标月和目标气压基本流场V分量的气候态
    
11. v_climatic = vwnd_mon.sel(time=vwnd_mon['time.month']==mon, level=p).mean(dim='time')
    
12. 
    
13. # 经纬度转为弧度制
    
14. lon_deg = Φ['lon'].copy()
    
15. lat_deg = Φ['lat'].copy()
    
16. lon_rad = np.deg2rad(lon_deg) * units('1')
    
17. lat_rad = np.deg2rad(lat_deg) * units('1')
    
18. 
    
19. # 科氏参数
    
20. f = mpcalc.coriolis_parameter(Φ['lat'])
    
21. # 目标月和目标气压基本流场的气候态
    
22. wind_climatic = mpcalc.wind_speed(u_climatic, v_climatic)
    
23. 
    
24. cosφ = np.cos(lat_rad)
    
25. 
    
26. # 位势的纬向偏差
    
27. Φ_prime = Φ - Φ_climatic.mean(dim='lon')
    
28. # 将需要对弧度制经纬度求偏导的量的坐标都换成弧度制经纬度
    
29. Φ_prime = Φ_prime.assign_coords({'lon': lon_rad, 'lat': lat_rad})
    
30. f = f.assign_coords({'lat': lat_rad})
    
31. cosφ = cosφ.assign_coords({'lat': lat_rad})
    
32. u_climatic = u_climatic.assign_coords({'lon': lon_rad, 'lat': lat_rad})
    
33. v_climatic = v_climatic.assign_coords({'lon': lon_rad, 'lat': lat_rad})
    
34. wind_climatic = wind_climatic.assign_coords({'lon': lon_rad, 'lat': lat_rad})
    
35. # 准地转流函数相对于气候场的扰动
    
36. Ψ_prime = Φ_prime / f
    
37. 
    
38. # 一顿偏导猛如虎
    
39. dΨ_prime_dλ = Ψ_prime.differentiate('lon')
    
40. dΨ_prime_dφ = Ψ_prime.differentiate('lat')
    
41. ddΨ_prime_ddλ = dΨ_prime_dλ.differentiate('lon')
    
42. ddΨ_prime_ddφ = dΨ_prime_dφ.differentiate('lat')
    
43. ddΨ_prime_dλφ = dΨ_prime_dλ.differentiate('lat')
    
44. 
    
45. # T-N波作用通量的水平分量公共部分
    
46. temp1 = p / p0 * cosφ / (2 * wind_climatic * earth_avg_radius**2)
    
47. temp2 = dΨ_prime_dλ * dΨ_prime_dφ - Ψ_prime * ddΨ_prime_dλφ
    
48. 
    
49. # T-N波作用通量的水平分量
    
50. fx = temp1 * (u_climatic / cosφ**2 * (dΨ_prime_dλ**2 - Ψ_prime * ddΨ_prime_ddλ) + v_climatic / cosφ * temp2)
    
51. fy = temp1 * (u_climatic / cosφ * temp2 + v_climatic * (dΨ_prime_dφ**2 - Ψ_prime * ddΨ_prime_ddφ))
    
52. 
    
53. # 把弧度制经纬度再换成角度制便于画图
    
54. lon = lon_deg.values
    
55. lon[lon>180] -= 360
    
56. fx = fx.assign_coords({'lon': lon, 'lat': lat_deg}).sortby(['lon', 'lat'])
    
57. fy = fy.assign_coords({'lon': lon, 'lat': lat_deg}).sortby(['lon', 'lat'])

复制代码

查看一下fx可以看到，计算出来的单位确实是m2/s2

简单画个图



如有错误或者疏漏，欢迎指正。

                               
登录/注册后可看大图

crabdudu

清清楚楚干净利落 像你学习

来自苹果APP客户端

Fernando

很好的教程，给你一个大赞

little卓忘密码

这个好像不是原创吧？

墨家大宝

little卓忘密码 发表于 2021-5-29 16:13
这个好像不是原创吧？

保证原创

任性小王子

不错不错，受教了

来自苹果APP客户端