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

Chwithin · 发表于 2020-3-19 11:53:22

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由 Chwithin 于 2020-8-2 22:19 编辑

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, 下载次数: 2111)

灰色节能君 · 发表于 2020-3-21 12:42:12

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

MeteoInfo · 发表于 2024-6-11 09:36:32

haochizhenhaoch 发表于 2024-6-11 09:06
那太好了，我现在也那么改的。ε和sor_index你怎么取的？尤其是ε？

具体可以参考此贴：https://www.bilibili.com/read/cv35290342/

MeteoInfo · 发表于 2024-6-9 17:39:11

本帖最后由 MeteoInfo 于 2024-6-9 17:44 编辑

haochizhenhaoch 发表于 2024-6-4 21:08
楼主你好，我最近正需要辐散风场来分析RWS，也用了沙的matlab程序。正巧看到你主贴里也说了对沙程序中迭代 ...

同意，dy2在前，dx2在后。根据本贴和沙贴在MeteoInfoLab中实现了类似功能，修改dx2, dy2顺序后没有南北拉长到边界的现象了。

freekiller · 发表于 2020-3-19 13:14:01

赞一个

小其其格 · 发表于 2020-3-19 14:28:48

手动点赞！！！

olislae · 发表于 2020-3-24 23:57:08

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

sealsu · 发表于 2020-3-31 14:48:21

试试看，我需要

wangzhaoying · 发表于 2020-4-1 09:23:03

谢谢楼主分享

Chwithin · 发表于 2020-4-2 20:27:35

本帖最后由 Chwithin 于 2020-8-2 22:05 编辑

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

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

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

Chwithin · 发表于 2020-4-2 20:30:00

本帖最后由 Chwithin 于 2020-4-2 20:37 编辑

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

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

sealsu · 发表于 2020-4-2 22:21:01

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

		自动登录	找回密码
密码			立即注册

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

登录后查看更多精彩内容~

评分