- 实习目的
熟悉迭代方法在气象中的应用,掌握稳定度指数的实际编程计算。
- 实习内容
1.已知2020年的6h再分析资料,要素场有温度、高度、相对湿度和水平风场。用迭代法求某一区域(10ºN-60ºN,60ºE-140ºE) 各层各日(各小组负责日期的第2日)20时的假相当位温,并给出在850和500hPa上相应的图及其此两层之差。
2.已知2020年的6h再分析资料,要素场有温度、高度、相对湿度和水平风场。求各层各日(各小组负责日期的第2日)20时区域平均(70ºE-100ºE) 、(100ºE-110ºE)和(110ºE-120ºE)在(10ºN-60ºN,1000-300hPa)面上的假相当位温和位温,并给出相应的图。
- 实习结果分析
3.1 计算的程序
(1)计算850hpa假相当位温、500hpa假相当位温、两层之差
%读取.nc资料数据%变量
T=ncread(['E:\天诊\','air.2020.nc'],'air');%温度T
H=ncread(['E:\天诊\','hgt.2020.nc'],'hgt');%高度H
RH=ncread(['E:\天诊\','rhum.2020.nc'],'rhum');%相对湿度RH
u=ncread(['E:\天诊\','uwnd.2020.nc'],'uwnd');%水平风场u
v=ncread(['E:\天诊\','vwnd.2020.nc'],'vwnd');%水平风场v
%(4维)索引变量-lon,lat,level,time
lon=ncread(['E:\天诊\','air.2020.nc'],'lon');
lat=ncread(['E:\天诊\','air.2020.nc'],'lat');
level=ncread(['E:\天诊\','air.2020.nc'],'level');
time=ncread(['E:\天诊\','air.2020.nc'],'time');
%10?N-60?N,100?E-140?E——925hPa气压上——6月9日——北京20时为——露点温度、凝结高度和凝结温度与时间
timez=(datenum(2020,6,9,12,0,0)-datenum('1800-01-01 00:00:0.0'))*24; %将世界12时转化为资料时间
%转换索引下标
time_z=find(time==timez);
level_z=find(level==850);
minlon=find(lon==100);
maxlon=find(lon==140);
minlat=find(lat==10);
maxlat=find(lat==60);
%取出所需资料(lon,lat,level,time)
TT=T(minlon:maxlon,maxlat:minlat,level_z,time_z);
HH=H(minlon:maxlon,maxlat:minlat,level_z,time_z);
RHRH=RH(minlon:maxlon,maxlat:minlat,level_z,time_z);
uu=u(minlon:maxlon,maxlat:minlat,level_z,time_z);
vv=v(minlon:maxlon,maxlat:minlat,level_z,time_z);
%-800hpa-转换索引
time_z=find(time==timez);
level_z=find(level==850);
minlon=find(lon==100);
maxlon=find(lon==140);
minlat=find(lat==10);
maxlat=find(lat==60);
%取出所需资料(lon,lat,level,time)
TT=T(minlon:maxlon,maxlat:minlat,level_z,time_z);
HH=H(minlon:maxlon,maxlat:minlat,level_z,time_z);
RHRH=RH(minlon:maxlon,maxlat:minlat,level_z,time_z);
%先计算露点——之前编程用的是摄氏度公式
RH=RHRH;tz=TT-273.15;e0=6.1078;a=17.269;b=35.86;
es=e0.*exp((a.*tz)./(273.16+tz-b));
ez=RH./100.*es;
for i=1:maxlon-minlon+1
for j=1:minlat-maxlat+1
while ez(i,j)-es(i,j)<0
tz(i,j)=tz(i,j)-0.01;
es(i,j)=e0.*exp((a.*tz(i,j))./(273.15+tz(i,j)-b));
end
td(i,j)=tz(i,j);
end
end
%再计算假相当位温
RH=RHRH;Tz=TT;Td=td+273.15;P0=1000.0;Pz=925;Z=HH;e0=6.1078;a=17.269;b=35.86;Cp=1004;g=9.8;L=2.5*10^6; %2.5×10^6J/kg
ez=e0.*exp((a.*(Td-273.15))./(Td-b));
seita_z=Tz.*(P0./Pz).^0.286;
qz=0.622.*ez./(Pz-0.378.*ez);
for i=1:maxlon-minlon+1
for j=1:minlat-maxlat+1
Pl=300;
Tl(i,j)=seita_z(i,j).*(Pl./P0).^0.286;
el(i,j)=e0.*exp((a*(Tl(i,j)-273.15))/(Tl(i,j)-b));
ql(i,j)=0.622.*el(i,j)/(Pl-0.378.*el(i,j));
while ql(i,j)-qz(i,j)<0
Pl=Pl+1;
Tl(i,j)=seita_z(i,j).*(Pl./P0).^0.286;
el(i,j)=e0.*exp((a*(Tl(i,j)-273.15))/(Tl(i,j)-b));
ql(i,j)=0.622.*el(i,j)/(Pl-0.378.*el(i,j));
end
seita_se_850(i,j)=seita_z(i,j).*exp(L*qz(i,j)/(Cp.*Tl(i,j)));
end
end
%500hpa
clear TT HH RHRH %清除之前数据
level_z=find(level==500);
RH=ncread(['E:\天诊\','rhum.2020.nc'],'rhum');%相对湿度RH
%取出所需资料(lon,lat,level,time)
TT=T(minlon:maxlon,maxlat:minlat,level_z,time_z);
HH=H(minlon:maxlon,maxlat:minlat,level_z,time_z);
RHRH=RH(minlon:maxlon,maxlat:minlat,level_z,time_z);
%计算露点——之前编程用的是摄氏度公式
RH=RHRH;tz=TT-273.15;e0=6.1078;a=17.269;b=35.86;
es=e0.*exp((a.*tz)./(273.16+tz-b));
ez=RH./100.*es;
for i=1:maxlon-minlon+1
for j=1:minlat-maxlat+1
while ez(i,j)-es(i,j)<0
tz(i,j)=tz(i,j)-0.01;
es(i,j)=e0.*exp((a.*tz(i,j))./(273.15+tz(i,j)-b));
end
td(i,j)=tz(i,j);
end
end
%计算假相当位温
RH=RHRH;Tz=TT;Td=td+273.15;P0=1000.0;Pz=925;Z=HH;e0=6.1078;a=17.269;b=35.86;Cp=1004;g=9.8;L=2.5*10^6; %2.5×10^6J/kg
ez=e0.*exp((a.*(Td-273.15))./(Td-b));
seita_z=Tz.*(P0./Pz).^0.286;
qz=0.622.*ez./(Pz-0.378.*ez);
for i=1:maxlon-minlon+1
for j=1:minlat-maxlat+1
Pl=300;
Tl(i,j)=seita_z(i,j).*(Pl./P0).^0.286;
el(i,j)=e0.*exp((a*(Tl(i,j)-273.15))/(Tl(i,j)-b));
ql(i,j)=0.622.*el(i,j)/(Pl-0.378.*el(i,j));
while ql(i,j)-qz(i,j)<0
Pl=Pl+1;
Tl(i,j)=seita_z(i,j).*(Pl./P0).^0.286;
el(i,j)=e0.*exp((a*(Tl(i,j)-273.15))/(Tl(i,j)-b));
ql(i,j)=0.622.*el(i,j)/(Pl-0.378.*el(i,j));
end
seita_se_500(i,j)=seita_z(i,j).*exp(L*qz(i,j)/(Cp.*Tl(i,j)));
end
end
%画图程序(1)
%两层之差
seita_se_500_850=seita_se_500-seita_se_850;
(2)计算1000-300hPa的北京20时的区域平均(70ºE-100ºE) 、(100ºE-110ºE)、(110ºE-120ºE)、(10ºN-60ºN)上的假相当位温和位温。
clearvars -except T H lon lat level time_z
rh=ncread(['E:\天诊\','rhum.2020.nc'],'rhum');%相对湿度RH
level_max=find(level==300);%水平下标
level_min=find(level==1000);
minlat=find(lat==10);%纬度下标
maxlat=find(lat==60);
lonlon=[70 100 110 120];
for i=1:4
lon_index(i)=find(lon==lonlon(i));%经度下标
end
%循环计算出图
for l=1:3
clearvars -except T H rh lon lat level time_z level_min level_max minlat maxlat lonlon lon_index l ave_se ave_z
minlon=lon_index(l);
maxlon=lon_index(l+1);
for k=level_min:level_max
level_z=k;
RHRH=rh(minlon:maxlon,maxlat:minlat,level_z,time_z);
TT=T(minlon:maxlon,maxlat:minlat,level_z,time_z);
HH=H(minlon:maxlon,maxlat:minlat,level_z,time_z);
%计算露点——之前编程用的是摄氏度公式
RH=RHRH;tz=TT-273.15;e0=6.1078;a=17.269;b=35.86;
es=e0.*exp((a.*tz)./(273.16+tz-b));
ez=RH./100.*es;
for i=1:maxlon-minlon+1
for j=1:minlat-maxlat+1
while ez(i,j)-es(i,j)<0
tz(i,j)=tz(i,j)-0.01;
es(i,j)=e0.*exp((a.*tz(i,j))./(273.15+tz(i,j)-b));
end
td(i,j)=tz(i,j);
end
end
%计算假相当位温
RH=RHRH;Tz=TT;Td=td+273.15;P0=1000.0;Pz=925;Z=HH;e0=6.1078;a=17.269;b=35.86;Cp=1004;g=9.8;L=2.5*10^6; %2.5×10^6J/kg
ez=e0.*exp((a.*(Td-273.15))./(Td-b));
seita_z=Tz.*(P0./Pz).^0.286;
qz=0.622.*ez./(Pz-0.378.*ez);
for i=1:maxlon-minlon+1
for j=1:minlat-maxlat+1
Pl=300;
Tl(i,j)=seita_z(i,j).*(Pl./P0).^0.286;
el(i,j)=e0.*exp((a*(Tl(i,j)-273.15))/(Tl(i,j)-b));
ql(i,j)=0.622.*el(i,j)/(Pl-0.378.*el(i,j));
while ql(i,j)-qz(i,j)<0
Pl=Pl+1;
Tl(i,j)=seita_z(i,j).*(Pl./P0).^0.286;
el(i,j)=e0.*exp((a*(Tl(i,j)-273.15))/(Tl(i,j)-b));
ql(i,j)=0.622.*el(i,j)/(Pl-0.378.*el(i,j));
end
seita_se(i,j)=seita_z(i,j).*exp(L*qz(i,j)/(Cp.*Tl(i,j)));
end
end
ave_se(:,k-level_min+1,l)=mean(seita_se,1);
ave_z(:,k-level_min+1,l)=mean(seita_z,1);
%%循环出图
% figure
% m_proj('lambert','lon',lonlon(l:l+1),'lat',[10 60]); hold on
% m_grid('linestyle','-','tickdir','out','linewidth',3); hold on
% m_contourf(lonlon(l):2.5:lonlon(l+1),60:-2.5:10,seita_se',15,'linestyle','none');
% hold on
% colormap(m_colmap('jet',10));%填充颜色搭配
% colorbar('Location','eastoutside');%颜色栏
% m_coast('color',[0 0 0]);%海岸线-黑色三原色
% title([num2str(l),'区域',num2str(level(k)),'hpa','假相当位温(K)']);
% hold on
% %保存当前窗口的图像
% saveas(gcf,[num2str(l),'区域',num2str(level(k)),'hpa','.jpg'])
% close gcf %不加会重叠
end%k
end%l
3.2 计算的绘图程序和描述文件
- 绘图程序
%画图程序1
figure
m_proj('lambert','lon',[100 140],'lat',[10 60]); hold on
m_grid('linestyle','-','tickdir','out','linewidth',3); hold on
m_contourf(100:2.5:140,60:-2.5:10,seita_se_500',15,'linestyle','none'); hold on
%m_contour(100:2.5:140,60:-2.5:10,seita_se','linestyle','none');
colormap(m_colmap('jet',10));%填充颜色搭配
colorbar('Location','eastoutside');%颜色栏
m_coast('color',[0 0 0]);%海岸线-黑色三原色
title('假相当位温(K)');
hold on
%画图程序2
figure
m_proj('lambert','lon',[100 140],'lat',[10 60]); hold on
m_grid('linestyle','-','tickdir','out','linewidth',3); hold on
m_contourf(100:2.5:140,60:-2.5:10,seita_se_500_850',15,'linestyle','none'); hold on
%m_contour(100:2.5:140,60:-2.5:10,seita_se','linestyle','none');
colormap(m_colmap('jet',10));%填充颜色搭配
colorbar('Location','eastoutside');%颜色栏
m_coast('color',[0 0 0]);%海岸线-黑色三原色
title('假相当位温(K)');
hold on
%第二问画图
%画图程序3
Y=cell(1,8);
for i=level_min:level_max
Y(1,i-level_min+1)=cellstr([num2str(level(i)),'hpa']);
end
c(1,:)=10:5:60;
X=cell(1,11);
for i=1:11
X(1,i)=cellstr([num2str(c(1,i)),'N']);
end
for i=1:3
figure
subplot(1,2,1)
contourf(ave_z(:,:,i)')
title([num2str(i),'区域假相当位温经向垂直剖面图'])
set(gca,'YTick',1:8)
set(gca,'YTicklabel',Y')
set(gca,'XTick',1:2:21)
set(gca,'XTicklabel',X)
set(gca,'XTickLabelRotation',45);
subplot(1,2,2)
contourf(ave_se(:,:,i)');
title([num2str(i),'区域位温经向垂直剖面图'])
set(gca,'YTick',1:8)
set(gca,'YTicklabel',Y')
set(gca,'XTick',1:2:21)
set(gca,'XTicklabel',X)
set(gca,'XTickLabelRotation',45);
saveas(gcf,[num2str(i),'区域经向垂直剖面图','.jpg'])
close gcf %不加会重叠
end
3.3 绘制图形、统计表格和相关分析
- 第一问运行结果如图所示
图1 850hpa(上)与500hpa(下)假相当位温分布示意图
图2 与500hpa与850hpa假相当位温之差分布示意图
1.第二问绘图结果如图所示:
区域1(70ºE-100ºE) 、区域2(100ºE-110ºE)和区域3(110ºE-120ºE)在10ºN-60ºN,1000-300hPa,假相当位温经向垂直剖面图和位温经向垂直剖面图如下:
- 经向垂直剖面图
各层假相当位温分布
2.第二问结果分析:
青藏高原及其东部部分地区,850hpa与500hpa假相当位温差值较大,当位温廓线出现随高度增加而减小的情况,则说明此时气块的内能非常大,温度较高,可能产生剧烈的上升运动,即对流不稳定。可分析出,不仅空气中含有不少的水汽,并且此区域上层气块不稳定,随着气块上升 有大量水汽凝结脱落。
3.4 存在的问题或遇到的问题或体会或小结
温度较其他人较低,查了几十遍发现数据取出时原始温度就较低。但未发现程序错误。不知道怎么回事~