[原创]matlab直接插值提取wrfout中变量的工具包（持续更新）

sfhstcn2

MATLAB现在可以直接处理wrfout的数据了！——wrfout_post工具包！

目前实现的功能：
垂直层次由eta坐标插值到p坐标或z坐标；
水平插值；
p坐标和z坐标下地形空白点外推；
一部分诊断量的计算。

建议用法：
将此压缩包解压后得到的文件夹放在matlab安装目录的toolbox下然后set path。

注意事项：
（1）目前程序还不完善，可能也存在一些问题，持续更新、改进、完善中。
（2）使用方法在包里的ReadMe里面。
（3）提供了一个testget.m来试运行（速度绝对不慢！）。
（4）matlab要2012 或者7.0以上的版本，就是能支持netcdf读取的。

使用说明：
用于提取wrfout数据的函数：GetVar，用法：
[ var1,var2 ] = GetVar( filename,varname,LON,LAT,LEV,T,extrap,inter_typ )
其中
filename为文件名；
varname为所提取的变量名（参考wrf说明书上关于wrfout各变量意义）；
LON和LAT为所要得到的格点数据的经纬度值，一般用meshgrid得到；
LEV为所要插值的气压层，p坐标下单位hPa，z坐标下单位为km，可以为1*n的数组。注意当使用外推（extrap=1或'T'或'TRUE'）时，LEV的值需从地表到高空排列，且length不宜超过模式eta的层次；
T为所提取的时次，1*n的数组，若设为'all'则提取所有时次；
extrap为外推选项，即对于插值后因为地形等原因产生的缺省值点进行外推插值使其获得数值，设为1或者'T'或'TRUE'则执行外推，设为其它则不外推。
inter_typ为插值模式，'p'为插值到p坐标下，'z'为插值到z坐标下。
var1为提取的主变量，var2为辅变量。当varname设为下表诊断量中所标示的值时，提取结果为标示结果；而当varname不为以上标注时，var1输出处理后所提取的变量，var2可输出所提取结果的时间，格式为1*nt（nt为提取的时间维长度）的struct变量。struct的格式为：
{ char year，month，day，hour，minute，second } time
每个子变量均为字符型，依次表示年月日时分秒。如：time(1).year，time(5).minute，均可以表示某时刻的年或分时间信息。


另外提供了几个常用诊断量的提取，直接将varname项设为如下即可：
'height'       位势高度（m）（仅用于p坐标插值）
'pressure'     气压（hPa）（仅用于z坐标插值）
'tk'           开氏温度（K）
'theta'        位温（K）
'thetase'      假相当位温（K）
'td'           露点温度（℃）
'rh'           相对湿度（%）
'slp'          海平面气压（hPa）
'dbz'          3-D相当反射率（dBZ）(var1为dbz，var2为maxdbz）
'maxdbz'       最大反射率（dBZ）
'wspd'         风速（m/s)（var1为wspd，var2为wdir）
'wdir'         风向（°）（同上）
'ws10'         10米风速（m/s)（var1为ws10，var2为wd10）
'wd10'         10米风向（°）（同上）
'mpv'          湿位涡（PVU）（var1为mpv1，var2为mpv2）


9月30日更新日志：
（1）将提取变量的函数改为GetVar；
（2）提供了两个供输出的变量（var1和var2）用于提取作为一对的诊断量（如dbz，wspd，mpv等），当varname不为以上标注时，var1输出处理后要提取的变量，var2可输出模拟结果的时间，格式为1*nt（nt为提取的时间维长度）的struct变量。
（3）实现对z坐标的垂直插值和外推。
（4）增加3个用于地球坐标下计算涡度散度和梯度的函数：m_cur,m_div,m_grad。
（5）增加可提取的诊端量：'pressure','slp','wspd','wdir','ws10','wd10','mpv'。具体意义见上。

sfhstcn2

帅气小葛 发表于 2018-4-10 07:27
楼主，我发现NCL里面计算散度的公式是dv(j,i) = (v(j+1,i)-v(j-1,i))/dy2(j)+ (u(j,i+1)-u(j,i-1))/dx2(j)- ...

谢谢指出！确实是有问题，应当基于球坐标系计算，那么旋度和梯度应该也有问题，正在加紧斧正！
错误修改：
在m_div函数的结尾加上这句：
dv = dv-v.*tan(LAT*pi/180)/R;
在m_curl函数的结尾加上这句：
cl = cl+u.*tan(LAT*pi/180)/R;

或者直接下载这个并替换掉原来的：

Paulin

楼主 Get_GridVar.m 里面 xlon(xlon<0)=360-xlon(xlon<0);    貌似写错了，应该是xlon(xlon<0)=360+xlon(xlon<0);

帅气小葛

楼主，我发现NCL里面计算散度的公式是dv(j,i) = (v(j+1,i)-v(j-1,i))/dy2(j)+ (u(j,i+1)-u(j,i-1))/dx2(j)- (v(j,i)/a)*tan(phi(j))，The last terms accounts for the convergence of the meridians on a sphere.你的公式里面貌似丢了最后一项

sfhstcn2

ppap 发表于 2017-4-13 16:37
不好意思，我没在wrfout的参数中找到tk这个参数

wrfout里面是没有tk，但可以通过wrfout里的另外一些诊断量直接给算出来，而这个工具包里面有计算这个的程序，因此直接设为‘tk’是可以提取温度的

sfhstcn2

lhaikun@163.com 发表于 2015-9-17 09:52
帮你顶顶 用户使用量还有待提升

拉去数值模式版如何

lhaikun@163.com

楼主不将压缩包内容截个图展示一下同时来个使用说明吗？使用说明最好用编辑中添加代码的方式直接写到帖子里，让大家方便看，要不提供的信息太少了

Aires

感谢LZ分享！～

zhenliangli

赞一个！楼主好样的！

zhenliangli

平常都是用Grads提取
这下长见识了！

kongfeng0824

谢谢分享               

rs06

这个好，顶楼主下

ooo

非常好的程序，下载下来用用。

lhaikun@163.com

帮你顶顶 用户使用量还有待提升