利用Python计算质量流函数
（MSF）表征哈德来环流
作者：石头人
邮箱：2205455617@qq.com
数据：ERA5数据（自行下载）
# 导入库import xarray as xrimport numpy as npimport matplotlib.pyplot as plt
# 读取数据f = xr.open_dataset("E:/tp/mon/Wind for HC.nc")print(f.variables.keys())lat = np.array(f.latitude.loc[60:-60])# The meridional-vertical MSF is an important metric ofthe HC (Oort and Yienger 1996).# The MSF can be easilyobtained via a mass-weighted vertical integral of the zonal mean meridional velocity.v = np.array(f.v.loc[f.time.dt.month.isin([7])].loc[:,:,60:-60,110:120])lev = np.array(f.level)
# 计算MSFva = v.mean(0).mean(2)coslat = (np.cos(lat*np.pi/180))g = 9.8 # 重力加速度a = 6400000 # 地球半径k1 = lat.shape[0]fic = np.zeros([10,k1])for j in range(0,10,1):for k in range(0,k1,1):f0 = np.trapz(va[0:j,k],lev[0:j],axis=0) # 积分项# ERA5的数据第一层是100hPa，有的数据第一层是1000hPa，积分的时候要注意fic[j,k] = 2*np.pi*a*f0*coslat[k]/gfi = fic/(10**10) # 除以10^10后在图上更方便标明数字
# 画图fig = plt.figure(figsize=(15,10), dpi=600)ax1 = plt.subplot(1, 1, 1)c1 = ax1.contourf(lat,lev,fi,levels=np.arange(-1.5,1.55,0.1), extend = 'both',zorder=0, cmap=plt.cm.RdBu_r)c2 = ax1.contour(lat,lev,fi,levels=np.arange(-1.5,1.55,0.1), extend = 'both',zorder=1, colors = 'k',)ax1.set_ylabel('High (hPa)',fontsize=18)ax1.set_xlabel('Latitude',fontsize=18)ax1.set_yscale('symlog') # 对数坐标ax1.set_yticks([200,300,500,850,1000])ax1.set_yticklabels(['200','300','500','850','1000'],fontsize=15)ax1.set_ylim(1000,200) # 保证低空到高空是1000 hPa到100 hPaax1.set_xlim(-60,60)ax1.set_xticks([-60,-45,-30,-15,0,15,30,45,60])ax1.set_xticklabels([r'60$^\degree$S',r'45$^\degree$S', r'30$^\degree$S',r'15$^\degree$S',r'0$^\degree$',r'15$^\degree$N',r'30$^\degree$N', r'45$^\degree$N',r'60$^\degree$N'],fontsize=15)position=fig.add_axes([0.1, 0.023, 0.8, 0.017])fig.colorbar(c1,cax=position,orientation='horizontal',format='%.1f',)plt.savefig('E:/HC.jpg', dpi=600, bbox_inches='tight')
声明：欢迎大家转载、转发本公众号的原创作品，可以联系小编（微信号lizhi2412963740）授权。
第八星系人造大气理论爱好者致力于记录自己科研的脚印，致力于与大家一起交流学习进步！
