- 积分
- 1272
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2012-3-11
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2015-6-20 18:27:08
|
显示全部楼层
- % function [ output_args ] = eqtdraw( input_args )
- %EQTDRAW Summary of this function goes here
- % Detailed explanation goes here
- % 计算假相当位温并画图
- clc;clear all
- ncdisp('D:\fnl\fnl_20150418_12_00.grib2.nc');
- prs=(ncread('D:\fnl\fnl_20150418_00_00.grib2.nc','lv_ISBL0'))/100
- t1=ncread('D:\fnl\fnl_20150418_00_00.grib2.nc','TMP_P0_L100_GLL0');
- rh1=ncread('D:\fnl\fnl_20150418_00_00.grib2.nc','RH_P0_L100_GLL0');
- lon=ncread('D:\fnl\fnl_20150418_12_00.grib2.nc','lon_0');
- lat=ncread('D:\fnl\fnl_20150418_12_00.grib2.nc','lat_0');
- S2=shaperead('D:\国家基础地理信息系统数据\国界与省界\bou2_4p.shp');
- bou1x=[S2(:).X];
- bou1y=[S2(:).Y];
- for j=1:26
- lev(1:360,1:181,j)=prs(j);
- end
- es1=(6.112.*exp(17.67.*(t1-273.15)./(t1-29.65)));
- q1=rh1.*(0.62197.*es1./(lev-es1))/100;
- e1=lev.*q1./(0.62197+q1)+1e-10;
- tlcl1=55.0+2840./(3.5.*log(t1)-log(e1)-4.805);
- theta1=t1.*(1000./lev).^(0.2854.*(1.0-0.28.*q1));
- eqt1=theta1.*exp(((3376./tlcl1)-2.54).*q1.*(1.0+0.81.*q1))-273.5;
- tt1=eqt1(103:109,61:65,:);
- gg1=mean(mean(tt1));
- gp1=squeeze(squeeze(gg1));
- plot(flipud(gp1(6:end)),flipud(prs(6:end)),'-k*')
- set(gca,'YDir','reverse')
- axis([50 120 100 1000])
- % print('-dpng','D:\7.png','-r300')
复制代码 |
|