MATLAB丨绘制速度势和辐散风图（含代码、注释、数据、文献等）

Chwithin

MATLAB版本：MATLAB R2019a
OpenGrADS版本：OpenGrADS2.2.1
NCL版本：NCL6.4.0

数据：
ERA-interim(Monthly Means of Daily Means)  2019年1月 1000hPa 1°× 1° 的u、v风数据
原网址下载：https://apps.ecmwf.int/datasets/data/interim-full-moda/levtype=pl/

基本原理：
已知u、v风场 → 计算得到散度D → 设定边界条件求解泊松方程得到速度势φ → 对速度势φ求梯度得到散度风

                               
登录/注册后可看大图

画图（代码见附件）：
2019年1月份1000hPa速度势和辐散风图

                               
登录/注册后可看大图

                               
登录/注册后可看大图

                               
登录/注册后可看大图

讨论：
本帖计算速度势和辐散风所使用的matlab代码来源于气象家园早前的帖子（沙颖凯-2014.3.2：http://bbs.06climate.com/forum.php?mod=viewthread&tid=20777，沙颖凯-2014.2.26：http://bbs.06climate.com/forum.php?mod=viewthread&tid=20665，感谢大佬的分享！），但由于原帖没有贴图，自己最近刚好也需要画速度势和辐散风，就尝试着画了一下，并分别使用NCL和GrADS内置的函数画了同样的图进行对比验证，可以发现3个软件画出来的图高低值中心基本是一致的，总体上NCL和GrADS画出来的比较接近，matlab画的图在边界上比较发散（应该跟数值差分求解泊松方程时所设的边界条件和迭代求解的次数有关）。
值得一提的是：
1.NCL和GrADS内置的函数只能根据全球的风场来计算全球范围的速度势，无法根据区域上的风场数据求解区域速度势；而本贴使用到的matlab函数则可以，同时，如果对边界条件进行更加合理的设定且对迭代过程进一步优化，则求解得到的势函数结果会更加准确，具体可以参考相关文献（见附件）。
2.感觉（沙颖凯-2014.3.2：http://bbs.06climate.com/forum.php?mod=viewthread&tid=20777）楼主给的chi_potential函数使用的迭代方法貌似将理查逊法和加速利布曼法混用了？根据文献《有限区域流函数和速度势的3种求解方法在分析台风Bilis暴雨增幅中的比较研究(任晨平等,2013)》，我重新修订了速度势、辐散风函数phi_Velocity_Potential（理查逊法）。

附件预览：

                               
登录/注册后可看大图

附件： phi_Chwithin.zip (38.12 MB, 下载次数: 1969)

2020-3-19 13:44 上传

点击文件名下载附件

其他参考推荐：
NCL入门教程：https://www.bilibili.com/video/av25379202（实不相瞒，本贴的NCL-phi图是自己第一次使用NCL画的图）
NCL: Divergent and Rotational Wind Components: http://www.ncl.ucar.edu/Applications/wind.shtml
GrADS Extension Library for Calculating Streamfunction/Velocity Potential: http://opengrads.org/doc/udxt/fish/fish.html
How to compute a streamfunction?: http://pordlabs.ucsd.edu/matlab/stream.htm（还没看懂作者的算法）

2020年8月2日更新！！！
2020年3月19日的贴子发出后，有人反馈辐散风的风场矢量画错了，经过摸索，现特此订正！感谢大家的反馈！
新的风场矢量画法参考了姚志雄博士在科学网上的博文（http://blog.sciencenet.cn/blog-419857-624284.html），感觉姚志雄博士的思路是对的。
将下面的新附件 plot_phi_new.m 替换上面压缩包附件中的 plot_phi.m 文件，即可重新出图。
关于风场的绘制，还可以参考我的另外一篇帖子：绘制引导气流回归场（风场矢量）并进行显著性检验

新附件： plot_phi_new.m (3.58 KB, 下载次数: 227)

2020-8-2 21:45 上传

点击文件名下载附件

运行结果：

                               
登录/注册后可看大图

小白分享，如有错漏，敬请指正，欢迎交流！

灰色节能君

楼主你好，我一直不知道怎么加风场箭头标签，看了你的matlab代码后发现了一些问题，你的第一张图风场标签设置的是20m/s，我把10°E,50°S的Uphi设为20，Vphi设为0，按理说10°E,50°S处的风场箭头长度应该和标签箭头长度一致，但结果相差比较大。 不知道楼主知不知道哪里出了问题，交流一下吧

freekiller

赞一个                  

小其其格

手动点赞！！！

olislae

楼主你好，我用你的代码计算速度势异常与辐散风异常，画出来的图高值中心辐散，低值中心辐合，正好反了，请问你知道是什么问题吗？请问求异常应该如何求？

sealsu

试试看，我需要

wangzhaoying

谢谢楼主分享

Chwithin

灰色节能君 发表于 2020-3-21 12:42
楼主你好，我一直不知道怎么加风场箭头标签，看了你的matlab代码后发现了一些问题，你的第一张图风场标签设 ...

2020年4月2日回复：
我试了一下，确实错了，但是我暂时还没想到解决办法

2020年8月2日回复：
经过摸索，已重新更正，请试试新的附件，看看有没有错误，谢谢你的反馈！

Chwithin

olislae 发表于 2020-3-24 23:57
楼主你好，我用你的代码计算速度势异常与辐散风异常，画出来的图高值中心辐散，低值中心辐合，正好反了，请 ...

是否有些地方的正负号没有调好，从而导致的相反呢？
异常的话，应该先算出风场的异常，然后直接用附件里的函数 phi_Velocity_Potential 算出速度势和辐散风，即为速度势异常和辐散风异常吧。

sealsu

楼主，我用ncep 日资料画 144X73 200hpa时 '1990/02/28' 会迭代出不来，一直循环。