- 积分
- 163
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-6-19
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 F117_ren_0 于 2014-8-5 11:36 编辑
首先感谢气象家园各位元老,学习和借鉴了很多有用的东西,具体名字就不再点出来了,再次表示感谢!读取FNL数据首先安装read_grib工具箱和M_map工具箱,具体安装方法请百度!
1,读取FNL数据
grib_struct=read_grib('F:\FNL\fnl_19990801_00_00_c',-1); %%%读取NCEP再分析资料
read_grib('F:\FNL\fnl_19990801_00_00_c','inv'); %生成变量目录表文档
grib_struct(13).gds;
Hight_500hPa=grib_struct(213).fltarray;
%Hight_1000_500hPa=Hight_1000hPa_01;
for i=1:91;
for j=1:181;
lon(i,j)=181-j;
lat(i,j)=91-i;
Hight(i,j)=Hight_500hPa((i-1)*360+j)/10.0;
end
end
这里是grib1格式,好像是07年以后是grib2格式,但读法一样,只是从排列发生变化,具体可以见%生成变量目录表文档%想读取的数据 Hight_500hPa=grib_struct(213).fltarray里,只要改动括号里的数字就ok了
2,计算散度
clc
clear all
%%%%%%%%%%%%%%%%%%%%读取NFL
grib_struct=read_grib('F:\FNL\fnl_20100311_00_00_c',-1); %%%读取NCEP再分析资料
read_grib('F:\FNL\fnl_20100311_00_00_c','inv'); %生成变量目录表文档
%grib_struct(13).gds;
U850=grib_struct(157).fltarray;
V850=grib_struct(158).fltarray;
for i=34:50;
for j=70:100;
lon(i,j)=100-j;
lat(i,j)=50-i;
u850(i,j)=U850((i-1)*360+j);
v850(i,j)=V850((i-1)*360+j);
end
end
u850=u850(34:end,70:end);
v850=v850(34:end,70:end);
%计算散度
R0=6371004; pi=3.1415926;
lat0=50; lon0=70;
for lati=50:-1:34
latm(50-lati+1)= distance(lat0,lon0,lati,lon0)/180*pi*R0;
end
for loni=70:1:100
lonm(loni-69)= distance(lat0,lon0,lat0,loni)/180*pi*R0;
end
[xx,yy]=meshgrid(lonm,latm);
div=divergence(xx,yy,u850,v850);
save D:\matt\div850 div xx yy
850hPa散度计算完成,这里的经纬度我选的是新疆的范围,如有需要可自行改动。
3、计算涡度平流和温度平流
首先读取出UV,涡度,温度
lon=70:100;
lat=34:50;
PI=3.1415926;
R0=6371004;
[dtx850,dty850]=gradient(w850,1,1);
dx=gradient(lon,1)*PI/180;
dy=gradient(lat,1)*PI/180;
wp850=-((u850.*dtx850)./(dx.'*cos(lat*PI/180))'+v850.*(dty850'*diag(1./dy))')/R0;
以上w850为850hPa的涡度,wp850为计算出来的850hPa的涡度平流
[dtx850,dty850]=gradient(t850,1,1);
tp850=-((u850.*dtx850)./(dx.'*cos(lat.*PI/180))'+v850.*(dty850'*diag(1./dy))')/R0;
以上t850为850hPa的温度,tp850为计算出来的850hPa的温度平流。
4、计算K指数
clc
clear all
%%%%%%%%%%%%%%%%%%K指数
% K=(T850-T500)+Td850-(T700-Td700)
%%%%%%%%%%%%%%%%%%
load D:\matt\Temp2013041712
load D:\matt\RH2013041712
E0=611.2;
%%%%t>0 a=7.5; b=237.3;
%%%%t<=0 a=9.5 b=265.5
for i=1:17;
for j=1:31;
if t850(i,j)>0;
a=7.5;
b=237.3;
else
a=9.5;
b=265.5;
end
qq1(i,j)=a*t850(i,j);
qq2(i,j)=b+t850(i,j);
end
end
qq=qq1./qq2;
Es=E0*(10.^qq);%%%%%%%%饱和水汽压
E=Es.*rh850.*0.01; %%%该相对湿度下水汽压
T1=b*log10(E./Es);
T2=a-log10(E./Es);
Td850=T1./T2;
save D:\matt\Td850 Td850
clear all
load D:\matt\Temp2013041712
load D:\matt\RH2013041712
E0=611.2;
%%%%t>0 a=7.5; b=237.3;
%%%%t<=0 a=9.5 b=265.5
for i=1:17;
for j=1:31;
if t700(i,j)>0;
a=7.5;
b=237.3;
else
a=9.5;
b=265.5;
end
qq1(i,j)=a*t700(i,j);
qq2(i,j)=b+t700(i,j);
end
end
qq=qq1./qq2;
Es=E0*(10.^qq);%%%%%%%%饱和水汽压
E=Es.*rh700.*0.01; %%%该相对湿度下水汽压
T1=b*log10(E./Es);
T2=a-log10(E./Es);
Td700=T1./T2;
save D:\matt\Td700 Td700
clear all
clc
load D:\matt\Temp2013041712
load D:\matt\Td850
load D:\matt\Td700
K=(t850-t500)+Td850-(t700-Td700);
save D:\matt\Kzhishu2013041712 K
首先要读出湿度和温度,我以上程序中load出来的就是存好的数据,后面就是计算了。
5、画图
load D:\matt\wdpll2010031312;
fnshp_P='D:\Program Files\toolbox\m_map\bou2_4m\bou2_4p.shp';% 省界
S = shaperead(fnshp_P, 'UseGeoCoords', true,...
'Selector',{@(NAME) strcmp(NAME,'新疆维吾尔自治区'), 'NAME'});
figure(1);
set (gca,'Position',[0.1,0.1,0.85,0.85], 'color','w')
lo=70:100;
la=34:50;
[lo,la]=meshgrid(lo,la);
wp850=wp850.*10.^10;
wp850(isnan(wp850)==1)= 0;
[x,y]=meshgrid(70:0.1:100,34:0.1:50);
KK=griddata(lo(:),la(:),wp850(:),x,y,'v4');
m_proj('Equidistant Cylindrical','lon',[70,100],'lat',[34,50]);
[cs,h]=m_contour(x,y,KK,'LineWidth',2);
set(h,'ShowText','on','LevelList',[-500 -350 -200 -100 -50 0 50 100 200 350 500]);
clabel(cs,h,'fontsize',5);
hold on;
shading flat;
m_grid('box','fancy');
%axis equal
hold on
m_plot(S.Lon,S.Lat,'-k','linewidth',2);
hold on
title('850涡度平流');
close(1)
其他的也一样,就不一一表示了。
|
评分
-
查看全部评分
|