- 实习目的
（利用迭代法计算温湿特征参量）熟悉迭代方法在气象中的应用，掌握露点和凝结抬升高度的实际编程计算。
- 实习内容
1.已知南京站某一时次的相对湿度（97%）、气温（18.4ºC）、气压（1003.5hPa）以及站点高度（35.2m），请利用迭代法求出该站的露点温度、凝结高度和凝结温度。
2.已知2020年的6h再分析资料，要素场有温度、高度、相对湿度和水平风场。用迭代法求某一区域(10ºN-60ºN,100ºE-140ºE)在925hPa气压上各日(各小组负责日期的第2日)20时的露点温度、凝结高度和凝结温度，并给出相应的图。
（要求: 根据实习内容和资料说明，编写计算你熟悉的语言程序并输出计算结果）
- 实习结果分析
3.1 计算的程序
（1）计算露点温度
RH=0.97;Z=35.2;tz=18.4;e0=6.1078;a=17.269;b=35.86;
es=e0*exp((a*tz)/(273.16+tz-b));
ez=RH*es;
while ez-es<0
tz=tz-0.01;
es=e0*exp((a*tz)/(273.15+tz-b));
if(ez-es>0)
break;%跳出循环
end
end
Td=tz+273.15
（2）计算凝结高度和凝结温度
RH=0.97;Tz=18.4+273.15;P0=1000;Pz=1003.5;Pl=300;Z=35.2;e0=6.1078;a=17.269;b=35.86;Cp=1004;g=10.0;tz=18.4;
ez=e0*exp((a*(Td-273.15))/(Td-b));
seita_z=Tz*(P0/Pz)^0.286;
qz=0.622*ez/(Pz-0.378*ez);
Tl=seita_z*(Pl/P0)^0.286;
el=e0*exp((a*(Tl-273.15))/(Tl-b));
ql=0.622*el/(Pz-0.378*el);
while ql-qz<0
Pl=Pl+1;
Tl=seita_z*(Pl/P0)^0.286;
el=e0*exp((a*(Tl-273.15))/(Tl-b));
ql=0.622*el/(Pl-0.378*el);
if(ql-qz>0)
break;%跳出循环
end
end
Zl=(g*Z+Cp*Tz-Cp*Tl)/g;
Zl,Tl
（3）运行结果如图所示：
3.2 计算的绘图程序和描述文件
- 读取nc文件
了解数据结构——再分析数据来自数据集'NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Model'，来源于（NCEP）美国国家环境预测中心。关键变量数据维度为4，对应索引变量为lon,lat,level,time（144x73x17x1464）。
气压垂直方向17层，单位为hpa，每层具体数据如上图所示。纬度[90 -90]，经度[0 357.5]分辨率为2.5度
数据以1800-1-1 00:00:0.0为时间原点，是matlab的系统时间，转化后发现时间是从2020年1月1日0点整刻开始，每6小时一次观测。为调用数据，使用时需要将2020年6月8日至11日转换为matlab日期序列值。
计算所需区域(10ºN-60ºN,100ºE-140ºE)在925hPa气压上，7班A组负责日期6月8-11日的6月9日（第二日）20时的露点温度、凝结高度和凝结温度与时间。
程序代码
%program WY_1_td
RH=0.97;Z=35.2;tz=18.4;e0=6.1078;a=17.269;b=35.86;
es=e0*exp((a*tz)/(273.16+tz-b));
ez=RH*es;
while ez-es<0
tz=tz-0.01;
es=e0*exp((a*tz)/(273.15+tz-b));
end
Td=tz+273.15;
% program WY_1_Zl_Tl
RH=0.97;Tz=18.4+273.15;P0=1000;Pz=1003.5;Pl=300;Z=35.2;e0=6.1078;a=17.269;b=35.86;Cp=1004;g=10.0;tz=18.4;
ez=e0*exp((a*(Td-273.15))/(Td-b));
seita_z=Tz*(P0/Pz)^0.286;
qz=0.622*ez/(Pz-0.378*ez);
Tl=seita_z*(Pl/P0)^0.286;
el=e0*exp((a*(Tl-273.15))/(Tl-b));
ql=0.622*el/(Pl-0.378*el);
while ql-qz<0
Pl=Pl+1;
Tl=seita_z*(Pl/P0)^0.286;
el=e0*exp((a*(Tl-273.15))/(Tl-b));
ql=0.622*el/(Pl-0.378*el);
end
Zl=(g*Z+Cp*Tz-Cp*Tl)/g;
Zl,Tl
%读取.nc资料数据
%路径
indir = 'E:\天诊\';
ncdisp([indir,'air.2020.nc']);
ncdisp([indir,'hgt.2020.nc']);
ncdisp([indir,'rhum.2020.nc']);
ncdisp([indir,'uwnd.2020.nc']);
ncdisp([indir,'vwnd.2020.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');
%save('first.mat','T','H','RH','u','v','lat','lon','level','time')
%观察数值
level
lon
%前三次观测时间
Time0=datestr(datenum('1800-01-01 00:00:0.0') + time(1)./24)
Time1=datestr(datenum('1800-01-01 00:00:0.0') + time(2)./24)
Time2=datestr(datenum('1800-01-01 00:00:0.0') + time(3)./24)
%将2020,6,8,0,0,0转换为系统时间
time0=(datenum(2020,6,8,0,0,0)-datenum('1800-01-01 00:00:0.0'))*24 %将资料原本时间转化为系统时间
%检查转化正确
Time=datestr(datenum('1800-01-01 00:00:0.0') + time0/24)
%%%——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);%2.5
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);
%计算露点——之前编程用的是摄氏度公式
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
%计算凝结高度和凝结温度Zl,Tl
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;
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
Zl(i,j)=(g.*Z(i,j)+Cp.*Tz(i,j)-Cp.*Tl(i,j))./g;
end
end
%Cubic双立方插值
[lat0,lon0]=meshgrid(60:-2.5:10,100:2.5:140);
[lat1,lon1]=meshgrid(60:-1:10,100:1:140);
TTT=interp2(lat0,lon0,TT,lat1,lon1,'cubic');%
[lat0,lon0]=meshgrid(60:-2.5:10,100:2.5:140);
[lat1,lon1]=meshgrid(60:-1:10,100:1:140);
TD=interp2(lat0,lon0,Td,lat1,lon1,'cubic');%
[lat0,lon0]=meshgrid(60:-2.5:10,100:2.5:140);
[lat1,lon1]=meshgrid(60:-1:10,100:1:140);
TL=interp2(lat0,lon0,Tl,lat1,lon1,'cubic');%
[lat0,lon0]=meshgrid(60:-2.5:10,100:2.5:140);
[lat1,lon1]=meshgrid(60:-1:10,100:1:140);
ZL=interp2(lat0,lon0,Zl,lat1,lon1,'cubic');%
[lat0,lon0]=meshgrid(60:-2.5:10,100:2.5:140);
[lat1,lon1]=meshgrid(60:-1:10,100:1:140);
HHH=interp2(lat0,lon0,HH,lat1,lon1,'cubic');%
[lat0,lon0]=meshgrid(60:-2.5:10,100:2.5:140);
[lat1,lon1]=meshgrid(60:-1:10,100:1:140);
RHRR=interp2(lat0,lon0,RHRH,lat1,lon1,'cubic');%
RHRR=RHRR';
HHH=HHH';
ZL=ZL';
TD=TD';
TL=TL';
TTT=TTT';
AA=RHRR(51:-1:1,:);
HHHH=HHH(51:-1:1,:);
TTTT=TTT(51:-1:1,:)-273.15;
ZZL=ZL(51:-1:1,:);
TTD=TD(51:-1:1,:)-273.15;
TTL=TL(51:-1:1,:)-273.15;
%作风场图
[lonlon,latlat]=meshgrid(100:2.5:140,60:-2.5:10);
uu=uu';
vv=vv';
m_proj('lambert','lon',[100 140],'lat',[10 60]);
[CS,CH]=m_etopo2('contourf',[-5000:500:0 250:250:3000],'edgecolor','none');
m_grid('linestyle','none','tickdir','out','linewidth',3);
colormap([ m_colmap('blues',80); m_colmap('gland',48)]);
brighten(.5);
ax=m_contfbar(1,[.5 .8],CS,CH);
title(ax,{'Level/m',''}); % Move up by inserting a blank line
m_quiver(lonlon,latlat,uu,vv,'k')
hold on3.3 绘制图形、统计表格和相关分析
观察风场发现，有东南气流携带太平洋的大量暖湿水汽到达我国东南沿海。在20N-30N，112E-120E的区域有风场辐合，利于水汽聚集。在120E-136E也有一个辐合中心。
等压面的形势是南高北低，观察温度场发现南侧高温，北侧低温，温度随纬度的增加而递减。温度高等压面的厚度就打，使得南侧等压面较高。
在120E-136E的辐合区域、在20N-30N，112E-120E的区域有风场辐合利于水汽聚集，对应地区相对湿度较高，空气接近饱和。
在40N附件区域相对湿度极低，空气难以达到饱和，露点温度相应较低。而在120E-136E的辐合区域、在20N-30N，112E-120E的区域相对湿度较高，空气接近饱和。且温度较高，故露点较高，接近温度。
同理，在40N附件区域相对湿度极低，空气难以达到饱和，凝结需要抬升的高度越高。而在120E-136E的辐合区域、在20N-30N，112E-120E的区域相对湿度较高，空气接近饱和。且个别高温高湿站点凝结高度接近站点高度。
3.4 存在的问题或遇到的问题或体会或小结
（1）编程语法的遗忘使得编程困难
（2）第一次实习报告计算抬升凝结高度有负值，但未检查出错误。几十次检查后犯下是循环时，计算的初始条件忘记赋初值，导致抬升凝结气压随纬度累减，减为负值。——后来设置在循环中设置初始参数，结果正确，故修改第一次报告。
4317
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
