grads脚本计算P坐标系下的涡度方程

冬日情愫+

最近接到一个任务是计算P坐标系下涡度方程中等号右边的各项。想了想貌似觉得grads实现起来应该挺方便了，于是开始着手写脚本。不过写之前也搜索了气象家园，看看有没有现成的脚本可供参考的，不得不说气象家园是个宝藏库，搜索到 http://bbs.06climate.com/forum.p ... 0%B6%C8%B7%BD%B3%CC 中楼主（桔梗花开）之前分享过的帖子，下载下来学习了一下，给了我很不错的启发，当然我也同意在这个帖子下面一位论坛作者（我爱问问题）指出对β项计算应该漏除以地球半径，基于这个帖子的基础上我自己写了计算P坐标系下的涡度方程脚本。

                               
登录/注册后可看大图

我这里使用的数据是ERA5的高空数据（2.nc），高空变量相对比较少，而脚本中用到的变量主要由U风，V风，W风 ，然后还有数据本身自带的相对涡度vo，散度d。当然假如使用的数据中没有直接的相对涡度和散度变量的话用u,v风也很容易的根据vo=hcurl(u,v) ，d=hdivg(u,v) 计算得到相应的涡度和散度，在脚本中，由于涡度方程各项都量级都比较小，所以等号右侧5项都乘以了10的8次方。
指出几点计算涡度方程中的注意事项吧，如果有不对的地方，请回帖批评指正：1、脚本中用的是时间和层次循环批量处理的方法，我下载的数据从1000hPa到200hPa，使用脚本计算的时候提醒注意的是，由于在垂直方向使用的是中央差分的方法，所以最低层次和最高层次是不能包括在循环中的，我这边输出的是常用的几个气压层数据，最底层输出的是925hPa；
iz=1
while(iz<=5)
if(iz=1) ;'set lev 500'; zz=500hPa; endif
if(iz=2) ;'set lev 600'; zz=600hPa; endif
if(iz=3) ;'set lev 700'; zz=700hPa; endif
if(iz=4) ;'set lev 850'; zz=850hPa; endif
if(iz=5) ;'set lev 925'; zz=925hPa; endif
3、另一个需要提醒的是，如果需要用色斑输出，建议大家定义变量的范围稍微大一点，实际图像输出的范围稍微小一点，以免出现图像白边。在图片输出这块，我只是选择最简单的等值线默认输出，如果想要更美观的输出，可以根据自己需要设置。当然如果脚本中计算有什么不对的地方，敬请指出。
3、根据涡度方程可知，涡度方程左侧的单位应该是s**-2，方程右边的相对涡度平流和地转涡度平流项以及散度项中计算应该没有问题，主要注意的是除以的地球半径应该是6.37e6，即单位是m，和风场的单位一直。值得注意的是垂直速度的单位pa/s，所以 偏p的处理的时候注意也应该以pa为单位的，具体的在(lev(z+1)-lev(z-1))*100.0得到的是pa为单位。

vorticity_equation_calculate.gs (3.46 KB, 下载次数: 65)

2021-11-9 10:01 上传

点击文件名下载附件

mello_123

正需要，下载看看

myj0219

谢谢分享！！！！！！！！！！！！！！