天气学诊断分析I 实习报告(三)

  • 实习目的

熟悉迭代方法在气象中的应用,掌握稳定度指数的实际编程计算。

  • 实习内容

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. 绘图程序

%画图程序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. 第一问运行结果如图所示

 图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 存在的问题或遇到的问题或体会或小结

             温度较其他人较低,查了几十遍发现数据取出时原始温度就较低。但未发现程序错误。不知道怎么回事~

  • 0
    点赞
  • 34
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值