MATLAB计算水汽通量和水汽通量散度程序的讨论

kongfeng0824

本人在网上看到有MATLAB计算水汽通量的程序，特别感兴趣。
目前论坛和网上做水汽分析大都用NCL和GrADS，而MATLAB则几乎很少，这对于用MATLAB的同学来说，可能不太方便。
由于对本人不是大气科学背景，因此，可能对以下程序理解上有误，特把程序贴在这里，让大家帮忙看看。
核心想解决的问题是：
1.计算水汽通量WVF
2.计算水汽通量散度Div-WVF
3.计算从东西南北不同方向的水汽输入输出


程序如下：


% http://blog.sina.com.cn/s/blog_4cfb5a620101cesh.html
% 资料：hgt.mon.mean.nc；uwnd.mon.mean.nc，vwnd.mon.mean.nc；shum.mon.mean.nc(比湿)；
% 数据下载：http://www.esrl.noaa.gov/psd/dat ... lysis.pressure.html
%% -------------------------------------------------------------------------------------------------------------------------------
clc    %清屏
clear  %清内存变量
cd('D:\BaiduYunDownload');
filename='hgt.mon.mean.nc';
ncdisp(filename)
hgt1 = double(ncread('hgt.mon.mean.nc','hgt',[9 13 1 1],[65 37 8 793],[1 1 1 1]));
lon1  = double(ncread('hgt.mon.mean.nc','lon',9,64,1));
lat1  = double(ncread('hgt.mon.mean.nc','lat',13,37,1));
j=1;
for i=6:12:740
    hgt(:,:,:,[j,j+1,j+2])=squeeze(hgt1(:,:,:,[i,i+1,i+2]));
    j=j+1;
end
logc1=hgt<0;
hgt(logc1)=0.0;
sum11=62*3-sum(logc1);
sum1=sum(hgt);
[ii]=find(sum11==0);
sum11(ii)=1;
hgtm=squeeze(sum1./sum11);
clear file f hgt hgt1 logc1 sum11 sum1 ii


%% -------------------------------------------------------------------------------------------------------------------------------
filename1='uwnd.mon.mean.nc';
ncdisp(filename1)
uwnd1 = double(ncread('uwnd.mon.mean.nc','uwnd',[9 13 1 1],[65 37 8 793],[1 1 1 1]));
j=1;
for i=6:12:740
    uwnd(:,:,:,[j,j+1,j+2])=squeeze(uwnd1(:,:,:,[i,i+1,i+2]));
    j=j+1;
end
logc1=uwnd>200;
uwnd(logc1)=0.0;
sum11=62*3-sum(logc1);
sum1=sum(uwnd);
[ii]=find(sum11==0);
sum11(ii)=1;
uwm=squeeze(sum1./sum11);
clear file f uwnd uwnd1 logc1 sum11 sum1 ii j


%%  -------------------------------------------------------------------------------------------------------------------------------
filename2='vwnd.mon.mean.nc';
ncdisp(filename2)
vwnd1 = double(ncread('vwnd.mon.mean.nc','vwnd',[9 13 1 1],[65 37 8 793],[1 1 1 1]));
j=1;
for i=6:12:740
    vwnd(:,:,:,[j,j+1,j+2])=squeeze(vwnd1(:,:,:,[i,i+1,i+2]));
    j=j+1;
end
logc1=vwnd>200;
vwnd(logc1)=0.0;
sum11=62*3-sum(logc1);
sum1=sum(vwnd);
[ii]=find(sum11==0);
sum11(ii)=1;
vwm=squeeze(sum1./sum11);
clear file f vwnd vwnd1 logc1 sum11 sum1 ii j


%% -------------------------------------------------------------------------------------------------------------------------------
filename3='shum.mon.mean.nc';
ncdisp(filename3)
shum1 = double(ncread('shum.mon.mean.nc','shum',[9 13 1 1],[65 37 8 793],[1 1 1 1]));
j=1;
for i=6:12:740
    shum(:,:,:,[j,j+1,j+2])=squeeze(shum1(:,:,:,[i,i+1,i+2]));
    j=j+1;
end
logc1=shum>2000;
shum(logc1)=0.0;
sum11=62*3-sum(logc1);
sum1=sum(shum);
[ii]=find(sum11==0);
sum11(ii)=1;
shm=squeeze(sum1./sum11);
clear file f shum shum1 logc1 sum11 sum1 ii j


%% -------------------------------------------------------------------------------------------------------------------------------
g = 9.8;
qu = squeeze(sum((hgtm(:,1:7,:)-hgtm(:,2:8,:)).*(uwm(:,1:7,:)+uwm(:,2:8,:)).*(shm(:,1:7,:)+shm(:,2:8,:))/4,2));
qv = squeeze(sum((hgtm(:,1:7,:)-hgtm(:,2:8,:)).*(vwm(:,1:7,:)+vwm(:,2:8,:)).*(shm(:,1:7,:)+shm(:,2:8,:))/4,2));
qu=qu/g;
qv=qv/g;
zz=sqrt(qu.*qu + qv.*qv);
[x,y] = meshgrid(lon1,lat1);


load coast
lonlim=[70 160];
latlim=[10 60];
[Lat,Lon]=meshgrid(lat1(1:1:end),lon1(1:1:end));
%% Figure
figure('Color', 'w')
axesm('MapProjection', 'Eqdcylin', 'MapLatLimit', latlim, 'MapLonLimit', lonlim);
framem on; gridm on; plabel on; mlabel('South')
setm(gca, 'MLineLocation', 10, 'PLineLocation', 10)
setm(gca, 'MLabelLocation', 10, 'PLabelLocation', 10)
setm(gca, 'fontsize', 20, 'fontname','times new roman')
tightmap; axis off; hold on
linem(lat, long, 'Color', 'k', 'LineWidth', 1.5)
% quiver
Handle=quiverm(Lat, Lon, -qv',-qu', 1.0); % <--------------------------- size: 2.0
set(Handle, 'Color', 'b', 'LineWidth', 1.0);
% hold on;
% [cs,hh]=contourfm(Lat, Lon, zz',0:100:max(zz),'linestyle','none');
% % contourcmap('jet', 'Colorbar', 'on',  'Location', 'vertical', 'TitleString', 'Contour Intervals in 100','fontsize', 20);
% colormap(jet(20));
% caxis([0 1350]);
% colorbar('location','eastoutside','fontsize', 20);
%% quiver label
labelx=[72 73];
labely=[12.5 13]; % location of label % <------- Use the actual grid distance of the data
labelu=[500 0];
labelv=[0 0]; % set 25m/s as a standard
scale_auto=Get_Autoscale(Lon, Lat, -qv',-qu');
scale_label=Get_Autoscale(labelx, labely,  labelu, labelv);
scale_factor=scale_auto/scale_label;
Handle=quiverm(labely, labelx, labelv, labelu, 2.0*scale_factor); % <----- size: 2.0*scale_factor
set(Handle, 'Color', 'r', 'LineWidth', 1.5)
textm(labely(1)+1.5, labelx(2), [num2str(labelu(1)) ' m/s'],...
    'FontWeight', 'bold', 'Color', 'r', 'FontName', 'times new roman', 'FontSize', 20)
title('\fontsize{14}Water Vapor Flux', 'FontWeight', 'bold', 'FontName', 'Helvetica')
print(gcf,'-r300','-djpeg','NCEP_Water Vapor Flux');


来自群组: 北京师范大学地学-气象

咸鱼

为啥子我感觉那个水汽通量被弄得这么复杂，书上的公式不是简化了之后只有风速*比湿/g就可以了~

又是那隻貓

不是特别懂这种~

denghaijun

先赞一个！

kongfeng0824

咸鱼 发表于 2016-6-5 15:17
为啥子我感觉那个水汽通量被弄得这么复杂，书上的公式不是简化了之后只有风速*比湿/g就可以了~

你提的对。

AnneChen

赞，前阵子都找疯了。

kongfeng0824

AnneChen 发表于 2016-6-5 15:55
赞，前阵子都找疯了。

这是讨论帖子。不是最终的正确版本。

kongfeng0824

denghaijun 发表于 2016-6-4 23:16
先赞一个！

看看是否正确？

kongfeng0824

又是那隻貓 发表于 2016-6-4 22:53
不是特别懂这种~

哦哈哈。没事儿

wlzhongouc

我不是气象的啊，没有算过水汽通量，不过之前有收藏过一个，粘在这里了
* This is a script for displaying moisture convergence
* Written by Michael Maxwell
*   
* rh    = Relative Humidity in %
* t     = Temp at *set level in degrees Kelvin
* tc    = Temp in degrees C
* td    = Dewpoint at *set level in degrees C
* e     = Vapor pressure
* mixr  = Mixing ratio
* u     = U-wind in m/s
* v     = V-wind in m/s
* mconv = moisture convergence/divergence. convergence is positive and divergence is negative.
'reinit'
'open d:\12090.ctl'
'set lev 700'
'set lat 20 35'
'set lon 100 120'
'set t 5'
'set mpdset cnworld'
'set xlopts 1 4 0.15'
'set ylopts 1 4 0.15'
'set grads off'
'set timelab off'
'set grid off'
#'tc=(tmpprs-273.16)'
'td=tc-((14.55+0.114*tc)*(1-0.01*rh) + pow((2.5+0.007*tc)*(1-0.01*rh),3) + (15.9+0.37*tc)*pow((1-0.01*rh),14))'
'vapr=6.112*exp((17.67*td)/(td+243.5))'
'e=vapr*1.001+(lev-100)/900*0.0034'
'define mixr=0.62137*(e/(lev-e))*1000/9.8'
#通过此gs文件算出的水汽通量的单位是g/(cm.s),一般为 十几到几十
'define qx=u*mixr'
'define qy=v*mixr'
'define mconv=hdivg(qx,qy)*1e6'
'enable print d:\2lb-shuiqi900.gmf'
'set arrscl 1 500'
#'set lon 124'
#'set cmin 30'
'set gxout shaded'
'd mag(qx,qy)'
#'set cmin 30'
'set gxout contour'
'd mag(qx,qy)'
#'d theta'
#'d skip(u,8);skip(v,8)'
'print'
'disable print'