实习4  利用风场和高度场计算运动特征参量

• 实习目的

熟悉地转风涡度在气象中的应用，掌握涡度和散度的实际编程计算。

• 实习内容

已知2020年的6h再分析资料，要素场有温度、高度、相对湿度和水平风场。求各小组负责第二日20时10-60N，60-160E的区域在200hpa和 850hpa上的涡度和散度以及500地转风涡度，并给出相应的图。

• 实习要求:

根据实习内容和资料说明，编写计算你熟悉的语言程序(如Fortran、c和matlab等)，并输出计算结果和绘图(如NCL, Matlab和grads)。实习结果分析。

3.1 计算的程序

%读取.nc资料数据

indir = 'E:\天诊\';%路径

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');

%6月9日北京20时

timez=(datenum(2020,6,9,12,0,0)-datenum('1800-01-01 00:00:0.0'))*24; %将世界12时转化为资料时间

%转换索引下标

%-200hpa

time_z=find(time==timez);

level_z=find(level==200);

minlon=find(lon==60);

maxlon=find(lon==160);

minlat=find(lat==10);

maxlat=find(lat==60);

%取出所需资料(lon,lat,level,time)

uu=u(minlon:maxlon,maxlat:minlat,level_z,time_z);

vv=v(minlon:maxlon,maxlat:minlat,level_z,time_z);

latt=lat(maxlat:minlat)/180*pi;

uu=uu';

vv=vv';

%计算涡度

r=6371*10^3;

a=1/(2*r);

d=2.5/180*pi;%经纬度网格距转为弧度制

for i=2:size(uu,1)-1

    for j=2:size(uu,2)-1

            vor(i,j)=a*((vv(i,j+1)-vv(i,j-1))/cos(latt(i))/d -(uu(i+1,j)-uu(i-1,j))/d+2*uu(i,j)*tan(latt(i)));

    end

end

%计算散度

for i=2:size(uu,1)-1

    for j=2:size(uu,2)-1

              D(i,j)=a*((uu(i,j+1)-uu(i,j-1))/(cos(latt(i))*d) -(vv(i+1,j)-vv(i-1,j))/d-2*vv(i,j)*tan(latt(i)));

    end

end

%对站点数据进行Cubic双立方插值

X=60+2.5:160-2.5;

Y=10+2.5:60-2.5;

lonlon=double(lon(minlon+1:maxlon-1));

latlat=double(lat(minlat-1:-1:maxlat+1));

 

%850hpa

level_z=find(level==850);

uu=u(minlon:maxlon,maxlat:minlat,level_z,time_z);

vv=v(minlon:maxlon,maxlat:minlat,level_z,time_z);

uu=uu';

vv=vv';

%计算涡度

for i=2:size(uu,1)-1

    for j=2:size(uu,2)-1

            vor(i,j)=a*((vv(i,j+1)-vv(i,j-1))/cos(latt(i))/d -(uu(i+1,j)-uu(i-1,j))/d+2*uu(i,j)*tan(latt(i)));

    end

end

%计算散度

for i=2:size(uu,1)-1

    for j=2:size(uu,2)-1

              D(i,j)=a*((uu(i,j+1)-uu(i,j-1))/(cos(latt(i))*d) -(vv(i+1,j)-vv(i-1,j))/d-2*vv(i,j)*tan(latt(i)));

    end

End

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%500hpa

H=ncread(['E:\天诊\','hgt.2020.nc'],'hgt');%高度H

HH=H(minlon:maxlon,maxlat:minlat,level_z,time_z);

level_z=find(level==500);

uu=u(minlon:maxlon,maxlat:minlat,level_z,time_z);

vv=v(minlon:maxlon,maxlat:minlat,level_z,time_z);

uu=uu';

vv=vv';

HH=HH';

%计算地转风涡度

r=6371*10^3;

d=2.5/180*pi;%经纬度网格距转为弧度制

for i=2:size(uu,1)-1

    for j=2:size(uu,2)-1

      a=9.8/(2*(7.292*10^-5)*sin(latt(i))*r^2);

dz_vor(i,j)=a*((HH(i,j+1)+HH(i,j-1)-2*HH(i,j))/cos(latt(i))^2/d^2+(HH(i+1,j)+HH(i-1,j)-2*HH(i,j))/d^2-(H(i+1,j)-H(i-1,j))*tan(latt(i))/(2*d));

    end

end

3.2 计算的绘图程序和描述文件

%200涡度

[X1 Y1 Z1]=griddata(lonlon,latlat,double(vor(2:20,2:40)),X',Y,'cubic');%范围缩小一圈是起点是+

m_proj('miller','lon',[60 160],'lat',[10 60]);

m_grid('linestyle','none','box','fancy','tickdir','out');hold on

m_contourf(X1,Y1,Z1,30,'linestyle','none'); hold on

title('200hpa涡度分布图', 'Rotation', 0, 'FontSize', 16);%显示图标题

colormap(m_colmap('jet',10));%填充颜色搭配

colorbar('Location','eastoutside');%颜色栏

m_coast('color',[0 0 0]);%海岸线-黑色三原色

hold on;

[ulon,ulat]=meshgrid([60:2.5:160],[10:2.5:60]);

m_quiver(ulon,ulat,uu,vv,'k')%风场

hold on

%200散度

[X1 Y1 Z1]=griddata(lonlon,latlat,double(D(2:20,2:40)),X',Y,'cubic');%范围缩小一圈是起点是+

m_proj('miller','lon',[60 160],'lat',[10 60]);

m_grid('linestyle','none','box','fancy','tickdir','out');hold on

m_contourf(X1,Y1,Z1,30,'linestyle','none'); hold on

title('200hpa散度分布图', 'Rotation', 0, 'FontSize', 16);%显示图标题

colormap(m_colmap('jet',10));%填充颜色搭配

colorbar('Location','eastoutside');%颜色栏

m_coast('color',[0 0 0]);%海岸线-黑色三原色

hold on;

[ulon,ulat]=meshgrid([60:2.5:160],[10:2.5:60]);

m_quiver(ulon,ulat,uu,vv,'k')%风场

hold on

%850涡度

figure

[X1 Y1 Z1]=griddata(lonlon,latlat,double(vor(2:20,2:40)),X',Y,'cubic');%范围缩小一圈是起点是+

m_proj('miller','lon',[60 160],'lat',[10 60]);

m_grid('linestyle','none','box','fancy','tickdir','out');hold on

m_contourf(X1,Y1,Z1,30,'linestyle','none'); hold on

title('850hpa涡度分布图', 'Rotation', 0, 'FontSize', 16);%显示图标题

colormap(m_colmap('jet',10));%填充颜色搭配

colorbar('Location','eastoutside');%颜色栏

m_coast('color',[0 0 0]);%海岸线-黑色三原色

hold on;

[ulon,ulat]=meshgrid([60:2.5:160],[10:2.5:60]);

m_quiver(ulon,ulat,uu,vv,'k')%风场

hold on

%850散度

figure

[X1 Y1 Z1]=griddata(lonlon,latlat,double(D(2:20,2:40)),X',Y,'cubic');%范围缩小一圈是起点是+

m_proj('miller','lon',[60 160],'lat',[10 60]);

m_grid('linestyle','none','box','fancy','tickdir','out');hold on

m_contourf(X1,Y1,Z1,30,'linestyle','none'); hold on

title('850hpa散度分布图', 'Rotation', 0, 'FontSize', 16);%显示图标题

colormap(m_colmap('jet',10));%填充颜色搭配

colorbar('Location','eastoutside');%颜色栏

m_coast('color',[0 0 0]);%海岸线-黑色三原色

hold on;

[ulon,ulat]=meshgrid([60:2.5:160],[10:2.5:60]);

m_quiver(ulon,ulat,uu,vv,'k')%风场

hold on

%500地转风涡度

figure

[X1 Y1 Z1]=griddata(lonlon,latlat,double(dz_vor(2:20,2:40)),X',Y,'cubic');%范围缩小一圈是起点是+

m_proj('miller','lon',[60 160],'lat',[10 60]);

m_grid('linestyle','none','box','fancy','tickdir','out');hold on

m_contourf(X1,Y1,Z1,30,'linestyle','none'); hold on

title('500hpa地转风涡度分布图', 'Rotation', 0, 'FontSize', 16);%显示图标题

colormap(m_colmap('jet',10));%填充颜色搭配

colorbar('Location','eastoutside');%颜色栏

m_coast('color',[0 0 0]);%海岸线-黑色三原色

hold on;

[ulon,ulat]=meshgrid([60:2.5:160],[10:2.5:60]);

m_quiver(ulon,ulat,uu,vv,'k')%风场

hold on

3.3 绘制图形、统计表格和相关分析