请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 959|回复: 13

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

[复制链接]

新浪微博达人勋

发表于 2020-3-19 11:53:22 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 Chwithin 于 2020-3-20 17:10 编辑

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, 下载次数: 241)
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2020-3-21 12:42:12 | 显示全部楼层
楼主你好,我一直不知道怎么加风场箭头标签,看了你的matlab代码后发现了一些问题,你的第一张图风场标签设置的是20m/s,我把10°E,50°S的Uphi设为20,Vphi设为0,按理说10°E,50°S处的风场箭头长度应该和标签箭头长度一致,但结果相差比较大。 不知道楼主知不知道哪里出了问题,交流一下吧
{IJ2L[WWJ}_E7MU)34MW14U.png
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2020-3-19 13:14:01 | 显示全部楼层
赞一个                  
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-3-19 14:28:48 | 显示全部楼层
手动点赞!!!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2020-3-24 23:57:08 | 显示全部楼层
楼主你好,我用你的代码计算速度势异常与辐散风异常,画出来的图高值中心辐散,低值中心辐合,正好反了,请问你知道是什么问题吗?请问求异常应该如何求?
屏幕快照 2020-03-24 23.56.31.png
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-3-31 14:48:21 来自手机 | 显示全部楼层
试试看,我需要
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-4-1 09:23:03 | 显示全部楼层
谢谢楼主分享
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-4-2 20:27:35 | 显示全部楼层
灰色节能君 发表于 2020-3-21 12:42
楼主你好,我一直不知道怎么加风场箭头标签,看了你的matlab代码后发现了一些问题,你的第一张图风场标签设 ...

我试了一下,确实错了,但是我暂时还没想到解决办法
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-4-2 20:30:00 | 显示全部楼层
本帖最后由 Chwithin 于 2020-4-2 20:37 编辑
olislae 发表于 2020-3-24 23:57
楼主你好,我用你的代码计算速度势异常与辐散风异常,画出来的图高值中心辐散,低值中心辐合,正好反了,请 ...

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

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-4-2 22:21:01 | 显示全部楼层
楼主,我用ncep 日资料画 144X73 200hpa时 '1990/02/28' 会迭代出不来,一直循环。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表