专栏名称: 气象学家

【气象学家】公众号平台为您解读最新气象科研进展、分享气象实用编程技巧、追踪气象即时资讯。欢迎加入气象AI和Python交流群以及气象博士群！与5W+的专业人士一起交流互动！

气象编程 | 流函数与速度势

气象学家 · 公众号 · · 2024-06-12 02:43

正文

第一时间获取气象科研资讯

『气象学家』公众号交流群

加入

流函数和速度势一直广泛应用于大气、海洋、环境科学等领域。流函数对应旋转风，常用于大气环流演变、大洋环流模拟、大范围水汽输送分析等方面。速度势对应辐散风，代表大气中的地转偏差运动，密切联系大气垂直运动，进而产生水汽相变和潜热释放，在热带天气系统分析、中小尺度动力学、山谷风和湖陆风等局地环流研究和污染扩散等领域尤为重要。

全球流函数、速度势的计算由于不涉及边界条件有唯一解析解，而有限区域流函数、速度势求解问题，需要对有边界条件限制下的泊松方程求解，常用的方法有迭代法、谱展开法、变分法、调和余弦法等。参考气象家园论坛里沙颖凯用加速利布曼法计算流函数和速度势的MATLAB代码( http://bbs.06climate.com/forum.php?mod=viewthread&tid=20777) ，在MeteoInfo 3.8.11版本中加入了相关功能，考虑到计算速度，迭代计算的代码（liebermanIteration）用Java语言实现（https://github.com/meteoinfo/MeteoInfo/blob/54476160dd02766523325260ad53b40c627475e9/meteoinfo-math/src/main/java/org/meteoinfo/math/meteo/MeteoMath.java），其它部分用Jython语言实现，放在meteolib包中（https://github.com/meteoinfo/MeteoInfo/blob/54476160dd02766523325260ad53b40c627475e9/meteoinfo-lab/pylib/mipylib/meteolib/calc/kinematics.py），函数名分别为stream_function和velocity_potential，参数为经度、纬度、风场U分量、风场V分量、最大迭代次数（缺省为1e6）、最小残差值（缺省为1e-7）、超张弛系数（缺省为0.2）。

利用沙颖凯贴中提供的MATLAB数据测试了流函数：

fn = 'D:/Temp/matlab/Mar_21_2009_00UTM.mat'data = np.loadmat(fn)u = data['u'][::-1]v = data['v'][::-1]lon = data['longitude'][::-1]lat = data['latitude'][::-1]
psi, Upsi, Vpsi = meteolib.stream_function(lon, lat, u, v, 1000, 1e-5)
axesm()geoshow('country')    levs = arange(-40, 40.1, 5)contourf(lon, lat, psi*1e-6, levs, edgecolor='gray', extend='both', cmap='BlueRed')colorbar(shrink=0.8) qq = quiver(lon, lat, Upsi, Vpsi, color='darkgray',    size=10, overhang=1, headwidth=4, headlength=4, antialias=True)quiverkey(qq, 0.85, 0.12, 20, bbox={'edge':True, 'fill':True})title('Stream function')xlim(90, 160)ylim(0, 60)

上图结果和帖子中的图形类似，区别主要由于迭代算法代码中dx2、dy2顺序不同导致，由于数据二维数组维度分别为纬度（y）和经度（x），计算残差时第一维应该除以dy2，第二维除以dx2。

速度势测试利用了此贴（ http://bbs.06climate.com/forum.php?mod=viewthread&tid=94948&extra=page%3D1&page=1 ）中的数据：

fn = 'D:/Temp/nc/2019_mon_u_v.nc'f = addfile(fn)u = f['U'][0,::-1]v = f['V'][0,::-1]lon1 = u.dimvalue(1)lat1 = u.dimvalue(0)lon, lat = meshgrid(lon1, lat1)
phi, Uphi, Vphi = meteolib.velocity_potential(lon, lat, u, v)
axesm()geoshow('country', edgecolor='darkgray')    levs = arange(-8, 8.1, 2)contourf(lon1, lat1, phi * 1e-6, levs, extend='both', cmap='BlueRed')colorbar() ss = 10qq = quiver(lon[::ss,::ss], lat[::ss,::ss], Uphi[::ss,::ss], Vphi[::ss,::ss],     size=50, overhang=1, headwidth=4, headlength=4, antialias=True)quiverkey(qq, 0.85, 0.12, 5, bbox={'edge':True, 'fill':True})title('Velocity potential')

可以看出速度势的高低值中心和数值量级与NCL、GRADS的结果相似，也没有原MATLAB代码南北端色块拉长的边界的现象。

也利用NCL网站上的示例数据（https://www.ncl.ucar.edu/Applications/wind.shtml）进行了测试：

fn = 'D:/Temp/nc/uv300.nc'f = addfile(fn)u = f['U'][0]v = f['V'][0]lon1 = u.dimvalue(1)lat1 = u.dimvalue(0)lon, lat = meshgrid(lon1, lat1)
phi, Uphi, Vphi = meteolib.velocity_potential(lon, lat, u, v)
axesm()geoshow('country', edgecolor='darkgray')    levs = arange(-8, 8.1, 1)contourf(lon1, lat1, phi * 1e-6, levs, extend='both', cmap='amwg256')colorbar() ss = 2qq = quiver(lon[::ss,::ss], lat[::ss,::ss], Uphi[::ss,::ss], Vphi[::ss,::ss],     size=150, overhang=1, headwidth=4, headlength=4, antialias=True)quiverkey(qq, 0.85, 0.12, 3, bbox={'edge':True, 'fill':True})title('Velocity potential')